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The CESE development is driven by a belief that a solver should (i) enforce conservation 
laws in both space and time, and (ii) be built from a non-dissipative (i.e., neutrally stable) 
core scheme so that the numerical dissipation can be controlled effectively. To provide a solid 
foundation for a systematic CESE development of high order schemes, in this paper we describe 
a new 4th-order neutrally stable CESE solver of the advection equation du/dt + adu/dx = 
0. The space-time stencil of this two-level explicit scheme is formed by one point at the 
upper time level and three points at the lower time level. Because it is associated with three 
independent mesh variables Uj, ( U x )j , and (uxx)™ (the numerical analogues of u, du/dx , and 
d 2 u/d: r 2 , respectively) and four equations per mesh point, the new scheme is referred to as the 
a(3) scheme. As in the case of other similar CESE neutrally stable solvers, the a(3) scheme 
enforces conservation laws in space-time locally and globally, and it has the basic, forward 
marching, and backward marching forms. These forms are equivalent and satisfy a space-time 
inversion (STI) invariant property which is shared by the advection equation. Based on the 
concept of STI invariance, a set of algebraic relations is developed and used to prove that the 
a( 3) scheme must be neutrally stable when it is stable. Moreover it is proved rigorously that 
all three amplification factors of the a( 3) scheme are of unit magnitude for all phase angles 
if \v\ < 1/2 (v = aAt/Ax). This theoretical result is consistent with the numerical stability 
condition \v\ < 1/2. Through numerical experiments, it is established that the a(3) scheme 
generally is (i) 4th-order accurate for the mesh variables w" and (u x )"; and 2nd-order accurate 
for (Uxx)j- However, in some exceptional cases, the scheme can achieve perfect accuracy aside 
from round-off errors. 


1. Introduction 

The space-time conservation element and solution element (CESE) method is a high-resolution and 
genuinely multidimensional method for solving conservation laws [1-73]. Its nontraditional features include: 

(i) a unified treatment of space and time; (ii) the introduction of conservation elements (CEs) and solution 
elements (SEs) as the vehicles for enforcing space-time flux conservation; (iii) a novel time marching strategy 
that has a space-time staggered stencil at its core and, as such, fluxes at an interface can be evaluated 
without using any interpolation or extrapolation procedure (which, in turn, leads to the method’s ability 
to capture shocks without using Riemann solvers); (iv) the requirement that each scheme be built from a 
non-dissipative core scheme and, as a result, the numerical dissipation can be controlled effectively; and (v) 
the fact that mesh values of the physical dependent variables and their spatial derivatives are considered as 
independent marching variables to be solve for simultaneously. Note that CEs are non-overlapping space- 
time subdomains introduced such that (i) the computational domain can be filled by these subdomains; and 

(ii) flux conservation can be enforced over each of them and also over the union of any combination of them. 
On the other hand, SEs are space-time subdomains introduced such that (i) the boundary of each CE can 
be divided into several component parts with each of them belonging to a unique SE; and (ii) within a SE, 
any physical flux vector is approximated using simple smooth functions. In general, a CE does not coincide 
with a SE. 

Without using flux-splitting or other special techniques, since its inception in 1991 [1] the unstructured- 
mesh compatible CESE method has been used to obtain numerous accurate ID, 2D and 3D steady and 
unsteady flow solutions with Mach numbers ranging from 0.0028 to 10 [51]. The physical phenomena 
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modeled include traveling and interacting shocks, acoustic waves, vortex shedding, viscous flows, detonation 
waves, cavitation, flows in fluid film bearings, heat conduction with melting and/or freezing, electrodynamics, 
MHD vortex, hydraulic jump, crystal growth, and chromatographic problems [3-73]. In particular, the rather 
unique capability of the CESE method to resolve both strong shocks and small disturbances (e.g., acoustic 
waves) simultaneously [13,15,16] makes it an effective tool for attacking computational aeroacoustics (CAA) 
problems. Note that the fact that second-order CESE schemes can solve CAA problems accurately is an 
exception to the commonly-held belief that a second-order scheme is not adequate for solving CAA problems. 
Also note that, while numerical dissipation is needed for shock capturing, it may also result in annihilation of 
small disturbances. Thus a solver that can handle both strong shocks and small disturbances simultaneously 
must be able to overcome this difficulty. 

In spite of its nontraditional features and potent capabilities, the core ideas of the CESE method are 
simple. In fact, all of its key features are the inescapable results of an honest pursuit driven by these 
simple ideas. The first and foremost is the belief that the method must be solid in physics. As such, in 
the CESE development, conservation laws are enforced locally and globally in their natural space-time unity 
forms for ID, 2D and 3D cases. Moreover, because direct physical interaction generally occurs only among 
the immediate neighbors, use of the simplest stencil also becomes a CESE requirement. Obviously, this 
requirement is also very helpful in simplifying boundary-condition implementation. 

The second idea emerges from the realization that stability and accuracy are two competing issues in 
time-accurate computations, i.e., too much numerical dissipation would degrade accuracy while too little 
of it will cause instability. In other words, to meet both accuracy and stability requirements, computation 
must be performed away from the edge (“cliff”) of instability but not too far from it. This represents a real 
dilemma in numerical method development. As an example, schemes with high-order accuracy generally has 
high accuracy and low numerical dissipation. However, it is susceptible to instability. In fact, in dealing with 
complicated real-world problems, stability of these schemes often is difficult to maintain without resorting 
to ad hoc treatments. To confront this issue head-on, in CESE development, it is required that a solver 
be built from a non-dissipative (i.e., neutrally stable) core scheme. By definition, computations involving 
a neutrally stable scheme are performed right on the edge of instability and therefore the numerical results 
generated are non-dissipative. As such numerical dissipation can be controlled effectively if the deviation of 
a solver from its non-dissipative core scheme can be adjusted using some built-in parameters. Note that the 
above idea also plays an essential role in the recent successful development of a family of Courant number 
insensitive schemes [59,61,64,65,67]. 

Other CESE ideas are: (i) the flux at an interface be evaluated in a simple and consistent manner; 
(ii) genuinely multidimensional schemes be built as simple, consistent and straightforward extensions of 
ID schemes; (iii) triangular and tetrahedral meshes be used in 2D and 3D cases, respectively, so that the 
method is compatible to the simplest unstructured meshes and thus can be used to solve problems with 
complex geometries; and (iv) logical structures and approximation techniques used be as simple as possible, 
and special techniques that has only limited applicability and may cause undesirable side effects be avoided. 
Fortunately for the CESE development, as it turns out, the realization of the above lesser ideas (i)-(iv) 
follows effortlessly from that of the first two core ideas. 

The first model equation considered in the CESE development is the simple convection equation 


du du 
~dt +a &c = 


( 1 . 1 ) 


where the advection speed a / 0 is a constant. Let Xi = x, and X 2 = t be considered as the coordinates 
of a two-dimensional Euclidean space En- Then, because Eq. (1.1) can be expressed as V • h = 0 with 

h = ( au,u ), Gauss’ divergence theorem in the space-time E 2 implies that Eq. (1.1) is the differential form 
of the integral conservation law 

</ h ■ ds = 0 (1.2) 

JS(v ) 

As depicted in Fig. 1, here (i) S(V) is the boundary of an arbitrary space-time region V in E 2 , and (ii) 
ds = da n with da and n, respectively, being the area and the unit outward normal vector of a surface element 
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on S'(U). Note that: (i) because h ■ ds is the space-time flux of h leaving the region V through the surface 
element ds, Eq. (1.2) simply states that the total space-time flux of h leaving V through S(V) vanishes; 
(ii) in E 2 , da is the length of a line segment on the simple closed curve 5(V); and (iii) all mathematical 
operations can be carried out as though E 2 were an ordinary two-dimensional Euclidean space. 

It is well known that a solution to Eq. (1.1) represents non-dissipative data propagation along its 
characteristic lines defined by dx/dt = a. Moreover, Eq. (1.1) is invariant under space-time inversion (STI), 
i.e., it transforms back to itself if x and t are replaced by — x and —t, respectively. (In physics, STI invariance 
generally is referred to as PT invariance where P denotes a mirror-image or spatial-reflection operation while 
T denotes a time-reversal operation). Thus a solution to Eq. (1.1) possesses the following properties: (i) it 
is completely determined by the data specified at an initial time level; (ii) its value at a space-time point 
has a finite domain of dependence (a point) at the initial time level; and (iii) the space-time inversion image 
of a solution to Eq. (1.1) is also a solution and vice versa. As such, in the initial CESE development, the 
focus is on the construction of an ideal core solver of Eq. (1.1) that enforces the conservation law Eq. (1.2) 
and also possesses all other properties of Eq. (1.1), i.e., it is a two-level, explicit, non-dissipative, and STI 
invariant solver. An in-depth account of this development and the resulting “a” scheme is given in [71]. As 
it turns out, the 2nd-order accurate a scheme (i) has a space-time stencil formed by one mesh point at the 
upper time level and two mesh points at the lower time level; and (ii) it is neutrally stable if v 2 < 1 where 
v = a At / ax. Also, at each space-time mesh point ( j,n ), the a scheme is associated with two independent 
mesh variables it" and ( u x )" (the numerical analogues of u and du/dx, respectively) and two equations. 

Until recently, with one exception (a three-level and 3rcl-order accurate scheme reported on p. 80 of [1]), 
all CESE solvers of Eq. (1.1) are two-level and 2nd-order accurate extensions of the a scheme. To initiate 
a systematic CESE development of high-order schemes, in this paper we describe a new 4th-order accurate, 
conservation-law enforcing, and neutrally stable CESE solver of Eq. (1.1). As will be shown, the space-time 
stencil of this two-level explicit scheme is formed by one point at the upper time level and three points at 
the lower time level. Because it is associated with three independent mesh variables it", (u x )" and («„)" 
(the numerical analogues of u, du/dx, and d 2 u/dx 2 , respectively) and three equations at each mesh point, 
hereafter the new scheme is referred to as the a(3) scheme. 

2. The a(3) scheme 

To proceed, consider the set fi of space-time mesh points ( j , n ) (marked by dots and crosses in Fig. 2(a)) 


where 

0 = f {(j,n)\j,n = 

0, ±1, ±2, ±3, . . .} 

(2.1) 

We have 

fi = fi 

1 U H 2 

(2.2) 

where Hi and H 2 are two disjoint sets defined by 



Hi = f {(j,n)\j,n = 0,±1,±2,±3,. 

. . , and (j + n ) is an odd integer} 

(2.3) 

^2 = f {(j, n)\j, n = 0,±1,±2, ±3, . 

. . , and (J + n) is an even integer} 

(2.4) 


In Fig. 2(a), the mesh points G Hi are marked by dots while those £ fl- 2 are marked by crosses. Hereafter 
O 2 is referred to as the complement set of Hi and vice versa. Obviously each of Hi and H 2 represents a set 
of space-time staggered mesh points. 

Each (j, n) € H is associated with (i) a solution element (SE), denoted by SE (j,n) (see Fig. 2(b) where 
(j, n) € Hi is assumed), and (ii) two conservation elements (CEs), denoted by CE ~{j,n) and CE + (j, n) (see 
Figs. 2(c) and 2(d) where (j, n) G Hi is assumed), respectively. Each SE is the interior of a space-time 
region that includes a horizontal line segment, a vertical line segment, and their immediate neighborhood. 
On the other hand, each CE is a rectangular space-time region. Hereafter, (i) SEs or CEs associated with 
mesh points € Hi (G H 2 ) may be referred to simply as SEs or CEs associated with Hi (H 2 ). 

As a preliminary for the following development, note that (see Figs. 2(a)-(d)): 
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(a) Two CEs which are associated with two mesh points, one £ fix while another £ fi 2 may occupy the 
same space-time region. As an example, (i) CE_(j, n) and CE + (j — 1 ,n) occupy the same space-time 
region; and (ii) (j, n) G fix (j — l,n) £ fi 2 . Hereafter the symbol “<t=>” is used as a shorthand for the 
statement “if and only if” . 

(b) A pair of diagonally opposite vertices of a CE both belong to the same set fix or fi 2 while another pair 
both belong to the complement set. As an example, points A and C belong to fix while points B and 
D belong to fi 2 . 

(c) The CEs associated with each of Ox and fi 2 by themselves are nonoverlapping and can fill the space-time 
£- 2 - 

(d) Among the line segments forming the boundary of the same space-time region occupied by both 
CE_(j, n) and CE+(j — l,n), (i) AB and AD C SE (j, n); (ii) CB and CD C SE (j — 1 ,n— 1); (iii) BA 
and BC C SE(j — l,n); and (iv) DA and DC C SE(j, n— 1). Because AB and BA represent the same 
line segment, one can see that any line segment on this boundary is a subset of two SEs with one of 
them being associated with Ox and another associated with 0 2 . Hereafter, this ambiguity is removed 
by the following SE designation rule: any line segment designated as a boundary of a CE associated 
with Ox (fi 2 ) is designated as a subset of a SE associated with Ox (0 2 ). As an example, if AB, AD, 
CB, and CD are designated as boundaries of CE_(j, n), then because points A and C belong to Ox, 
the above rule implies that: (i) both AB and AD are designated as subsets of SE(j, n); and (ii) both 
CB and CD are designated as subsets of SE(j — 1, n — 1). On the other hand, if BA, BC, DA, and DC 
are designated as boundaries of CE + (j — 1 ,n), then: (i) both BA and BC are designated as subsets of 
SE(j — 1, n)\ and (ii) both DA and DC are designated as subsets of SE(j, n — 1). 

Let ( x,t ) £ SE(j, n). Then Eqs. (1.1) and (1.2) will be simulated numerically assuming that u(x,t) and 
h{x,t), respectively, are approximated by 


u*{x, t ; j, n) 


and 


^u^ + iu^ix-x^ + iu^it-n+^u^ix-xjf + iu^ix-x^t-n + ^iuu^it-n 2 

(2.5) 

h*(x, t ; j, n) d = ( au*(x,t-j,n ), u*(x,t-,j,n )) (2.6) 


Note that: (i) it™, ( u x )", («*)!•', ( u xx )", (it^t)™, and (w t t)" are constants in SE(j, n), and the numerical 
analogues of the values of u, du/dx, du/dt, d 2 u/dx 2 , d 2 u/dxdt, and d 2 u/dt 2 at the mesh point ( j,n ), 
respectively; (ii) ( Xj,t n ) are the coordinates of the mesh point (j, n) where Xj = j ax and t n = nAt ; (iii) 
u*(x,t;j,n) represents a 2nd-order Taylor’s approximation of u; and (iv) Eq. (2.6) is the numerical analogy 
of the definition h = ( au,u ). 

For any (j, n) £ fi, let u = u*(x,t ; j,n) satisfy Eq. (1.1) for all ( x,t ) £ SE(j, n). Then one has 


(u t )j = -ci{u x )j, (u xt )] = -a(u xx )j, and (u tt )" = a 2 {u xx )™ , ( j,n)£fl (2.7) 

Substituting Eq. (2.7) into Eq. (2.5), one has 


u*(x,t\j,n ) = u] + ( u x )” [(x - Xj) - a(t- t n )] + ]r(u xx )j [(x - Xj) -a(t- t n )] 2 , 


(j, n) £ fi (2.8) 


i.e., u", ( u x )", and ( u xx )" are the only independent mesh variables associated with (j, n). 

With the above preliminaries, next we derive the flux conservation relations that underline the a(3) 
scheme. 

2.1. Flux conservation relations 

Let the flux of h* conserve over all CEs, i.e., 


h* ■ ds = 0, 


S(CE-(J») 


{j, n ) € 


(2.9) 
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Because (i) with respect to CE_(j, n), the outward unit normal vectors n at AB , AD, CD, and CB are 
(0, 1), (1,0), (0, —1), and (—1,0), respectively; and (ii) with respect to CE + (j, n), the vectors n at AF, AD, 
ED, and EF are (0,1), (—1,0), (0,-1), and (1,0), respectively, by using (i) the definitions given following 
Eq. (1.2), (ii) the above SE designation rule, and (iii) Eqs. (2.6) and (2.8), it can be shown that Eqs. (2.9) 
and (2.10) are equivalent to 


(1 + v) u — (1 — v)u s + 
and 

(1 — v) [w + (1 + v)llx + 


2(1 -u + u 2 ) 

3 


Uxx = (1 + v) u + (1 - v)u s + 

- i 


2(1 — v + is 2 ) 


2(1 + v + v 2 ) 


Uxx = (1 - v) U - (1 + v)u s + 


2(1 + v + v 2 ) 


respectively. Here: (i) v = a At, /ax is the Courant number; (ii) 

(Ux)]= and («xx)? d = ^f-(u xx )] (2.13) 

and (iii) to simplify notation, in the above and hereafter we adopt a convention that can be explained using 
an expression on the left side of Eq. (2.11) as an example, i.e. , 


n— 1 

5 

(j, n) e ft 

3~ 1 

(2.11) 

n— 1 

? 

(j, n) G fi 

3+ 1 

(2.12) 


(2.13) 


, \ 2(1 + v + I/ 2 ) 

U + (1 + V)Ux H Uxx 


= u] + (1 + u)(u s )] + 2(1 + " + " 2) (Uxx)? 


At this juncture, note that: 

(a) Because 

du AX du d 2 U (Ax) 2 8 2 U _ def x 

dx 2 dx an dx 2 4 dx 2 1 2 ax/2 

the normalized parameters (u-xY/ and («xx)", respectively, can be interpreted as the numerical analogues 
of the values at (j, n) of the first and second derivatives of u with respect to the normalized coordinate 
x. 

(b) By definition, points B and D depicted in Fig. 2(c) do not belong to either SE(j, n) or SE(j — 1, n — 1). 
This fact, however, does not pose a problem for flux evaluation over S(CE_(j,n)) because the values 
of h* at isolated points do not contribute to the flux of h* over a finite line segment. Similarly, the fact 
that points D and F depicted in Fig. 2(d) do not belong to SE(j, n) and SE(j + 1, n — 1) does not pose 
a problem for flux evaluation over S(CE+(J,n)). 

(c) According to the SE designation rule, each line segment such as AB depicted in Fig. 2(c) can be 
assigned with two different fluxes of h* , one is associated with il i (hereafter referred to as the fli-flux) 
and another associated with iY (hereafter referred to as the f^-fhix). As such, among those local 
conservation relations Eqs. (2.9) and (2.10), those associated with ( j,n ) € fii are completely decoupled 
from those associated with (j, n) € D 2 . Because Eqs. (2.9) and (2.10) are equivalent to Eqs. (2.11) and 
(2.12), respectively, it follows that each of the two systems of equations defined by Eqs. (2.11) and (2.12) 
is formed by two decoupled subsystems, one is associated with il i while another associated with il 2 - 

(d) Moreover, because (i) the vector h* at any interface separating two neighboring CEs associated with the 
same set Q-i (O 2 ) is evaluated using the information from the same SE, and (ii) the unit outward normal 
vector on the surface element pointing outward from one of these two neighboring CEs is exactly the 
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negative of that pointing outward from another CE, one concludes that the flux leaving one of these CEs 
through the interface is the negative of that leaving another CE through the same interface. Due to this 
interface flux cancelation and the fact that the CEs associated with each of Oi and tt -2 by themselves 
are nonoverlapping and can fill the space-time E 2 , the local conservation relations Eqs. (2.9) and (2.10) 
associated with (j,ri) £ fli (( j,n ) G S+ 2 ) lead to a global conservation relation, i.e., the total fii- (H 2 -) 
flux of h* leaving the boundary of any space-time region that is the union of any combination of CEs 
associated with the same set fli (ft 2 ) vanishes. 

Let 1 — v 2 0, i.e. 1 + v ^ 0 and 1 — v ^ 0. Then Eqs. (2.11) and (2.12) reduce to 


u — (1 — v)u s + 


2(1- v + v 2 ) 


2\ 1" 


. + (1 — v)lix + 


2(1 — v + v 2 ) 


, (j, n) £ Cl (2.14) 


3- 1 


and 


u + (1 + u)ux + 


2(1 + u + u 2 ) 


2\ 


u (IT v)ux + 


2(1 + u + v 2 ) 


, (j, n) £ ft (2.15) 


l+i 


respectively. Obviously, each of the two systems of equations defined by Eqs. (2.14) and (2.15) is also formed 
by two decoupled subsystems. Moreover, each component equation in Eq. (2.14) represents a stronger 
condition than the corresponding equation in Eq. (2.11) in the sense that the former implies the latter 
for any given v while the latter implies the former only if an extra condition (i.e., v ^ —1 for this case) 
is imposed. Similarly, each component equation in Eq. (2.15) represents a stronger condition than the 
corresponding equation in Eq. (2.12). These stronger conditions will be used in the construction of the a(3) 
scheme. 

As a preliminary to a later development, next we will take a side tour and introduce the concept of 
invariance under space-time inversion. 

2.2. Invariance under space-time inversion 

Let u = u(x,t) be a solution to Eq. (1.1) in the domain —00 < x, t < + 00 , i.e., 



du(x, t) du(x , t ) n 

dt +a dx ’ 

—00 < x, t < +00 

(2.16) 

Let 

x' = f — x and 

-to 

111 

(2.17) 

and 

u( x, t ) = f u(- 

-x, -t) 

(2.18) 

Then (i) Eq. (2.16) ++ 

du(x',t') | du(x',t') n 

dt 1 +a dx' ~ ’ 

—00 < x' , t' < +00 

(2.19) 

and (ii) 

d d 

— — = — — and 
dt' dt 

d d 

dx' dx 

(2.20) 

Thus Eq. (2.16) 

dii(x, t) du(x, t) 

at + “ ax = 0 ' 

— OO < x, t < +OO 

(2.21) 

In other words, if u = 

u(x,t) is a solution to Eq. (1.1), so must be u = u(x,t) and vice versa. 

Because the 

one-to-one mapping 

(x,t) <-> (-x,-t), 

— OO < X, t < +OO 

(2.22) 
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represents a space-time inversion (STI) operation, hereafter (i) a pair of functions such as u and u will be 
referred to as the STI images of each other; and (ii) a partial differential equation (PDE) such as Eq. (1.1) 
is said to be STI invariant if the STI image of a solution is also a solution and vice versa. 

Next let 


,(M) 


(*,*) = 


def d k+i u(x,t) 


and u^ k ’ e \x,t) = 


def d k+l u{x,t) 


dx k dt ( v ’ ' dx k dt t ’ 

Then, with the aid of the chain rule, Eqs. (2.17), (2.18), and (2.23) imply that 


— oo < x, t < + 00 ; k,£ = 0,1,2,... 

(2.23) 


t ) = d k+( u{-x,-t) = k+e d k+e u{x\t') 

dx k dt e K ’ dx ,k dt ,e -oo<x,t <+oo; M = 0,l,2,... (2.24) 

= (-1 ) fc + V M V , f ') = (-1 )* + *u ( W(-x,-t) 


i.e., 


u (k ’V(x,t) 


u ( k ^) (—x, —t) if (, k + £) is even 
- U (M) (_ ;B) - t ) if (k + £) is odd 


(2.25) 


According to Eq. (2.23), t/ 0,0 ) = u and u ( 0, °) = u. Thus Eq. (2.18) is a special case of Eq. (2.24) with 
k = £ = 0. 

In the following, the concept of STI invariance will be introduced for the a(3) scheme. As a preliminary, 
note that: (i) 

(j, n) <-> {-j, -n) (2.26) 


is the numerical analogue of the STI mapping Eq. (2.22); and (ii) w", (w x )", (u t )", (u xx )™, ( u xt )", and (u tt )" 
are the numerical analogues of the values of u, du/d x, du/dt, d 2 u/dx 2 , d 2 u/dxdt , and d 2 u/dt 2 , at the mesh 
point ( j,n ), respectively. Thus, motivated by Eq. (2.25), the one-to-one mapping 


( 0 ?~-( 0 =?; (« 0 ? ~ -(«*)=? 

(w xx )j (w xx )_"; (w x t )™ (u x t)_"; (' Uttfj 

is taken as the numerical analogue of the one-to-one mapping 

u^ k ’ e \x, t) <-> ii/ k ' l \x, t), — oo < £, t < + 00 ; fc,£=0,l,2,3 


(2.27) 


(2.28) 


For the independent mesh variables, by using Eq. (2.13), Eq. (2.27) reduces to 


(%)" ^ -(%)-" | , 

{Uxx)/j J \(u xx )_ n 


(j, n) G ft 


(2.29) 


Eq. (2.29) can be expressed as 


gtj,n) <-> U q(—j, —n), 


(j, n) G ft 


(2.30) 


where 


and 


9(j, «) = f ( («*)" 

(%x)j 


(j, n) G ft 


,10 0 

t/ d = | 0 -1 0 

0 0 1 


(2.31) 


(2.32) 
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The matrix U is unitary. In fact it is a real matrix with 


U = U - 1 (2.33) 

Hereafter (i) M -1 denotes the inverse of any nonsingular square matrix M; (ii) for each ( j,n ), Uq(—j, —n) 
is referred to as the STI image of q(j, n); and (iii) the set formed by Uq(—j, — n), ( j , n) € 0 is also referred to 
as the image of the set formed by q(j, n), ( j , n) G O. According to Eq. (2.33), q(j, n) = UUq(—(—j), — (— n)). 
Thus q(j,n ) is the STI image of Uq(—j,—n ) as an individual (j,n) or as the set defined over O. In the 
following, we will show that by itself each of the four subsystems of equations associated with Eqs. (2.14) 
and (2.15) is STI invariant, i.e., the subsystem maps onto an equivalent subsystem under the mapping 
Eq. (2.29). 

As an example, consider the subsystem of equations formed by the component equations associated with 
fii in Eq. (2.14). Let it be denoted as Eq. (2.14a). Under the mapping Eq. (2.29), Eq. (2.14a) maps onto 

-(n- 1) 

, (j,n)efl i (2.34) 

— O’— i) 

At this juncture, note that, in addition to changing the sign of each u, x, mapping Eq. (2.29) requires that the 
upper and lower indices j, n, j — 1, and n— 1 in Eq. (2.14a) be replaced by their negatives, respectively. This 
is different from simply replacing the symbols j and n everywhere with —j and —n, respectively. Moreover, 
to simplify argument, hereafter system B is referred to as the STI image of system A if A maps onto B 
under the mapping Eq. (2.29), e.g., the subsystem Eq. (2.34) is the STI image of Eq. (2.14a). Let 

j* = f 1 — j and n* ^ 1 — n, (j, n) G Hi (2.35) 

Then, by using the fact that ( j * + n*) + ( j + n) = 2 and therefore ( j*,n *) G Hi <+> (j, n) G Oi, Eq. (2.34) 
can be cast into the form 

n* — 1 

, (f,n*)Gf h (2.36) 

i*-i 

By comparing Eqs. (2.14a) and (2.36), one can see that the subsystem Eq. (2.14a) is identical to its STI 
image Eq. (2.34) (which is identical to Eq. (2.36)). Thus, under the mapping Eq. (2.29), Eq. (2.14a) maps 
onto itself, i.e., the subsystem Eq. (2.14a) is STI invariant. QED. 

The STI invariance of another three subsystems associated with Eqs. (2.14) and (2.15) can be established 
in a similar manner. As such the system formed by all component equations in each of Eqs. (2.14) and (2.15) 
is STI invariant. 

The three mesh variables at any (j, n) G H are linked to those at (j — l,n — 1) and (j + l,n — 1) by 
two component equations in Eqs. (2.14) and (2.15), respectively. In order that the three mesh variables at 
(j, n ) can be determined in terms of those mesh variables at the (n — l)th time level, in the next subsection 
we introduce an extra STI invariant condition that links the mesh variables at (j, n) with those at the mesh 
point ( j , n — l). 

2.3. A family of STI invariant solvers 

Consider the following system of equations: 

[u + aux + fdu^ = [u - au s +/3w 5 *]" _1 , (j,n) G Cl (2.37) 

where a and (3 are parameters independent of (j, n). By definition, (j, n) G fii (fl 2 ) <t=> ( j,n — 1) G fU (Hi). 
Thus, unlike Eqs. (2.14) and (2.15), the mesh variables associated with Oi are linked to those associated with 
fi 2 through Eq. (2.37). However, as will be shown, like a subsystem associated with Eq. (2.14) or Eq. (2.15), 
the system of equations Eq. (2.37) is STI invariant for any pair of a and (3. 


\ 2(1 — v + u z ) 

U — (1 — v)Ux H Uxx 


, 2(1 — v + zA) 

U + (1 - V)U S H Uxx 


. , 2(1 — v + zA) 

U + (1 - v)Ux H Uxx 


, 2(1-1/ + lA) 

u - (1 - v)ux H U S x 



The STI image of the system Eq. (2.37) is 

[u - au x + Pu ss ] Ij =[u + aux + , (j, n) G 0 (2.38) 

Let 

jf def an j n ' j (j] 77.) g O (2.39) 

Then because ( j',n ') G f i <+> (j, n) G fh Eq. (2.38) can be cast into the form 

[m + aux + @Uxx\jr =[u- au s + fe]", _1 , (/ , n') G (2.40) 

By comparing Eqs. (2.37) and (2.40), one can see that the system Eq. (2.37) is identical to its STI image 
Eq. (2.38) (which is identical to Eq. (2.40)). Thus, under the mapping Eq. (2.29), Eq. (2.37) maps onto 
itself, i.e., the system Eq. (2.37) is STI invariant. QED. 

Because each of Eqs. (2.14) and (2.15) is STI invariant, one can see that, for any pair of a and /?, the 
system formed by Eqs. (2.14), (2.15), and (2.37) is STI invariant. 

Next, the three mesh variables at any ( j , n) G will be solved in terms of those at (j — 1 , n — 1) , (j,n— 1) 
and (j + l,n — 1) using Eqs. (2.14), (2.15), and (2.37). Let 

A d = Ul + av) -2/3 (2.41) 

O 

and assume a/ 0. Then it can be shown that Eqs. (2.14), (2.15), and (2.37) <+■ 
n 4 r n i n— l 

Uj — — pX OtUx 1 P'^xx\ j 


1 

+ - 
A 

1(^-7 


2ai r . 2(1 + v + u 2 ) i™- 1 

u (1 + v)u x + _ u xx 

3 J L 3 Jj+i 

(j, n) G ft 

(2.42) 

1 

+ - 
A 


?) (l + 1/) + 

2a] r 7, , 2(1 — is + v 2 ) ]"-i 

U + (1 V)U X + Uxx 

0 J L O i j- 1 




, ‘tv r n V 

\ u x )j = 7 ^ [u - au s + fjUxx J , 


1 [-2(1- v + v 2 ) 
+ A [ 3 

1 r 2(1 + v + v 2 ) 
A l 3 


f3 u - (1 + v)u s H ^ 

P u + (1 - v)ux + — — ^ 


2(1 + ^ + tA) -j"- 1 


2(1 -v + v 2 ) ]«-! 

O ^XX 

3 U - 1 


(«ii)i = [w - " 

J A J 


1 — v + a 


1 + v — a 


u—{ l + v) 
u + (1 - v) 


2(1 + is + v 2 ) l"- 1 

Ux T ^ ^xx 

3 Jj+1 

2(1 -v + is 2 ) l"- 1 

U X T ~ 

3 -1,7-1 


(j,n)Gfi (2.43) 


(j.n)Gfi (2.44) 


For any pair of a and /3 with A / 0, Eqs. (2.42)-(2.44) represent a solver for Eq. (1.1). In the next 
subsection, we pick out the pair of a and P with which the solver will have the smallest truncation error 
(i.e., the highest order of truncation error) for Eq. (2.42). 

2.4. A study of truncation error 

Because, at each ( j,n ), Eqs. (2.42)-(2.44) represent a system of three equations for three independent 
mesh variables, Eqs. (2.42)-(2.44) represent a numerical analogue of a system of three coupled partial differ- 
ential equations (PDEs) with three dependent variables. (Eq. (1.1) is one of these PDEs). As such, in the 
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following study, three different symbols u, v, and w will be used to denote the analytical versions of u”, and 
the non-normalized variables (u x )" and (M ra )", respectively. Specifically, let u(x,t ), v(x,t), and w(x,t ) be 
functions having all the derivatives needed. Thus one can define 


v(x,t) = v(x,t) 


du(x , t) 
dx 


and w(x,t) = w(x,t) 


d 2 u(x, t) 
dx 2 


Also, as an example, one can define 


\dx e dt m 


def d e+m U . 

= X Or , — (jAx.nAt) 
dx e dt m 


f, m = 0,1,2,... 


(2.45) 


(2.46) 


Next we will consider the “analytical” version of Eq. (2.42) which results from replacing (i) it”, (w®)™, 
and (uxx)™, respectively, with u”, v ", and tc", for each ( j,n ); and (ii) the index n with n+1 everywhere. 
By using Eq. (2.13) and the fact that (j. n + 1) G -O- (j, n) G fi, the analytical form can be expressed as 


def 1 f ~n+l 4 


__ +i __ 


aAx _ 
u — X - 


1 r / 2ai^ 

" A [v ~3 

1 r/ 2av 
" a K~3 
(j,n) G Q; a/0 


2a 

“ TJ L 

, 2 a 

(! + !/) + — 1 L 


/3(ax) 2 


(1 + I/) AX _ 

u v - 


(1 + v + iz 2 )ax 2 n 

r w 

6 Jj+i 


(2.47) 


(1 — u) AX _ (1 — V + v 2 )ax 2 


6 


-w 


J i-i 


= 0 


By applying Taylor’s formula, it can be shown that 


. def f (du du\ 4(a — v)duAx fd 2 u ?d 2 u\At 

(< A ={(a + “&) + ^^&If + (ae-“S?) 2- 


1 |"2ai^ 3 
aL 3 


+ P{l-v 2 ) 


2v(v — a) d 2 u (Aa’) 2 (d 3 u (At) 2 a d 2 v (ax) 3 

3a dx 2 At 


1 


2a 


3A ^ (1 + ^ + ^.^ jgx 4< 

rd 4 u 4<9 4 un (At) 3 1 r2ai^ 3 


V ctt 3 ' <9x 3 / 6 3 a dx 2 At 

dw (ax) 3 a(l + 4zA) — 3/3^ — 2i/ 3 <9 3 u (ax) 3 


V dt 4 


dx 4 ) 24 6a L 3 


+ P{l-v 2 ) 


9a dx 3 At 

d 3 v (ax) 4 (3 d 2 w (ax) 4 


1 T2A 4 „ 2 

— [— + (1 + 2^ 

1 r2a 


■>('- T) 


C>X 3 At 

d 4 u (ax) 4 


<9x 4 At 


6a dx 2 At 
O [(At) 4 ] 


+ — — (1 - V + u 2 ) + /3(1 - v) 0[(Ax) 5 ]/At + 0 [(ax) 4 ] + 0[At(Ax) 3 ] 


1 


A + v + v 2 ) - /3(1 + I/)J [o[(Ax) 5 ]/At + 0 [(ax) 4 ] + 0[At(AX’) 3 ] 

(j, n) G Q, ; a ^ 0 


dv ( AX ) 5 
<9x At 


(2.48) 


Note that (ei)" defined in Eq. (2.47) is normalized by the factor (1/a t) so that the lowest-order terms in 
the above Taylor’s expansion contain the leading term ( du/dt + adu/dx) which is independent of At and 
ax. Also, in Eq. (2.48) a term is denoted by 0[(At) fl (ax)^ 2 ] if there exists a constant C > 0 and two fixed 
integers t\ > 0 and £2 > 0 such that the absolute value of this term < C(At)^ 1 (ax )^ 2 for all sufficiently small 
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At and ax. Note that, in determining the order of magnitude of a term such as O [(aj:) 5 ] in Eq. (2.48), 
the parameters a and (3 are not assumed to be constants independent of At and Ax. In fact, to reduce the 
truncation error of the a(3) scheme, they will be chosen to be functions of v (see Eqs. (2.58)) and thus vary 
with the ratio At/ Ax. 

In the following, let u = u(x, t), v = v(x, t), and w = w(x, t) be a solution to the system of PDEs formed 
by Eq. (1.1) and 


i.e., 


dv, d 2 u 

v — — = 0 and w — ■— r = 0 
ox dx z 


du, du du d 2 u 

— + a— =0, u - — = 0, and w - — - = 0 
dt dx dx dx z 


(2.49) 


(2.50) 


In other words, here the scheme Eqs. (2.42)-(2.44) is considered as a solver of the system of PDEs Eqs. (1.1) 
and (2.49). Eqs. (2.45) and (2.50) imply that 


Qi+mfi d i+m w 

= 0 and — — r— — = 0 


dx e 9t n 

d 2 u 


dx e dt r 


t,m = 0,1,2,... 


2 1 d 2 u _ / d d \ f du du'. _ 
dt 2 a dx 2 V dt a dx ) \ dt + a dx , 


d 3 u 

W 


5 d 3 u 


d 2 


d 2 


dt. 2 Q dtdx 


, 

dx 2 


du du'. 

a + “5J |s0 


,1'., _ u 

dt . 4 a dx 4 


d 2 


d 2 


d d \ { du du'. 


v dt. 2 dx 2 / 

Note that the first equation in Eq. (2.50), and Eqs. (2.52)-(2.54) are all special cases of 


(2.51) 

(2.52) 

(2.53) 

(2.54) 


i-\-m 

dx e dt m 


^ + (-1 ) fe - V — 

dt k dx k 


= 0 , 


, ?tj. = 0,1,2,...; A: =1,2,3... 


(2.55) 


With the hint provided by Eqs. (2.52)-(2.54), Eqs. (2.55) can be proved using the first equation in Eq. (2.50) 
and elementary algebra. 

By using Eqs. (2.46) and (2.51)-(2.54), one can see that (ei)" reduces to 


(ei)? = 


4 (a — v) du ax 2v(v — a) d 2 u (ax) 2 a(l + 4^ 2 ) — 3/3^ — 2v 3 d 3 u (ax) 3 


3 a dx At 


3a dx 2 At. 


1 


2v 4 ^ „ o 2az/\i d 4 u (ax) 4 '' n 




2a, 


dx 4 At. 


9a 

-0[(At) 4 ] 


dx 3 At, 


1 
A 

1 

a L 3 
(j, n) e Cl ; a ^ 0 


A 1^ — (1 - V + v 2 ) + /3(1 - I/)J |^0[(AX’) 5 ]/At + 0[(ax) 4 ] + 0[At(AX’) 3 ] 
^(1 + V + V 2 ) - /3(1 + u)} |"o[(AX’) 5 ]/At + 0 [(ax) 4 ] + 0[At(AX’) 3 ] 


(2.56) 


By definition, the expression on the right side of Eq. (2.56) represents the truncation error of Eqs. (2.42) 
if the scheme Eqs. (2.42)-(2.44) are considered as a solver of the system of PDEs Eqs. (1.1) and (2.49). 
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Here the values of a and (3 will be chosen so that the truncation error will reach the highest order. From 
Eq. (2.56), one can see that the coefficients of the three lowest-order terms in the truncation error vanish if 

a — is = 0 and a(l + 4is 2 ) — 2>/3v — 2v 3 = 0 


For the case is ^ 0, Eq. (2.57) <+> 


a = is and (3 = 


1 + 2is 2 


(2.57) 


(2.58) 


Next the a(3) scheme will be defined as the special solver with a and (3 being chosen according to Eq. (2.58). 

2.5. The basic and forward marching forms of the a(3) scheme 
Assuming Eq. (2.58), Eqs. (2.37), (2.41)-(2.44) and (2.56) reduce to 


l + 2v 2 

u + isu s H ugg 


1 + 2 v 2 n " 1 

u - vug H u S g 


w" = 2 u — isug 


1 + 2z/ -in-1 \ v 

3 ugg — 


a = 2/3 

2(1 + IS + V 2 ) -\n-l 

U - (1 + is)ug H Ugg 

3 Jj+i 


l~v 


2(1 — is + is 2 ) i n 1 

u + (1 - is/ug H Ugg 

o 1 

1 — 2is 


(2.59) 

(2.60) 

(j,n)Gfi (2.61) 


( \n o T 1 + 2 V in—l 

( Ug)j =2is[u- isug H Ugg \ 


, , 2(1 + is + is 2 ) 

u - (1 + v)ug H Ugg 


- n— 1 

- 3 + 1 


1 + 2 v 


. 2(1 — Z/ + Z/ 2 ) -ini 

+ (1 V)Ux ~ W’xx 

3 


/ \77 o r 1 ~\~ 2z/ 2 -in—i 

\U j xx) j — o[u ISUx + “ U xx\j ~ 2 l 


/i \ 2(1 + z/ + z/ 2 ) 

U - (1 + z/)iz* H w** 


■ n— 1 

- i+i 


3 

2 L 


2(1 — is + ^ 2 ) i "-- 1 

u + (1 - i/)«s H Ugg 

o 


(j, n) G O 

(2.62) 

0»Gfi (2.63) 


and 


( ei )?= — ( — 


24 \ dx A ) , 


(ax) 4 


At 


+ 2a 2 At(Ax) 2 + a 4 (At) 3 + O [(At) 4 ] + O [(ax) 4 ] 


+ O [At ( ax) 3 ] + O [(At) 2 (Ax)^ 


0[(a.t) e 

At 


(j,n)€fl (2.64) 


Note that: (i) the forms of the last four terms in Eq. (2.64) have been simplified using the definition 
is = aAt/ Ax ; and (ii) the expression on the right side of Eq. (2.64) represents the truncation error of 
Eq. (2.61) if the scheme formed by Eqs. (2.61)-(2.63) is considered as a solver of the system of PDEs 
Eqs. (1.1) and (2.49). 

Next we convert Eq. (2.62) into its analytical form by replacing (i) w”, (u x )/j, and ( u xx ) r - , respectively, 
with m™, u™, and tu", for each ( j,n ); and (ii) the index n with n+ 1 everywhere. By using (i) Eq. (2.13), (ii) 
is = aAt/ Ax, and (iii) the fact that (j, n + 1) G ft <+> (j. n) £ ft, then after a normalization by the factor 1/2, 
the analytical form can be expressed as 


(eo) n = f -v n+1 - 
[2>3 2 3 (ax) 2 

. 2aAt,\ 

2ax V ax / L 

1 / 2aAt,\ 

2ax V ax J L 


2aAt 


— (l 
\ax V 


aAt, (ax) 2 + 2a 2 (At) 2 i n 

U -~ V+ 12 H, 

Ax + aAt, , (ax) 2 + aAt , ax + a 2 (At) 2 J\ n 

u v + w 

2 6 J j+ 1 


(j, n) G H (2.65) 


Ax — aAt, _ (ax) 2 — a At, ax + a 2 (At) 

u H v + - — 1 — —w 


Jf-i 


= 0 
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Similarly, the analytical form of Eq. (2.63) can be expressed as 


( def 1^ +1 + 


6 * 
1 

(as) 2 

1 

(as) 2 


(as) 2 


aAt. „ (as) 2 + 2a 2 (At) 2 _ 

U-—V+ w 




Ax + aAt. „ (Ax) 2 + aAtAX + a 2 (At) 2 „ 

u v H w 

2 6 

As — aAt _ (ax) 2 — aAt.Ax + a 2 (At) 2 _ 

u H v + - — — —w 


j + 1 


(j, n)esi 


- 0-1 


= 0 


By using Taylor’s formula and Eq. (2.45), Eqs. (2.65) and (2.66) imply that 


< 77 —1— 

r dv dv d 

n 1 

/ du 

du \ i 

At 

|_ 

d 2 v (At) 2 

h 

to 

<5> 

r 

. dt dx dx 

\~dt 

dx) - 

2 + 

dt 2 4 

dx 2 

dw [a 2 (Af) 2 - (ax) 2 ] 

d 

(d 2 u 2 d 2 u' 

\ (At) 2 

(i + ^ 2 ; 

dx 

6 

dx 

\ dt 2 

dx 2 . 

/ 4 

12 

0[(M) 

I 3 ] + O [(At) 2 Ax] 

+ o 

[a1(ax) 2 

]+o[ 

(AX) 3 ] 

U, n) 


Tax) 2 


and 


(2.66) 


(2.67) 


(dv 

ra 2 

(du 

du\ 


.dx 2 

K dt. 

+ a d^J 


dt 


„ dw d 2 
2a — — T 3a 


dx 


dx 2 


At, 

~6 


dx 2 7 dt, 2 


dx 2 J 12 \dx 3 dx 2 7 6 


c> 2 /9 2 a 2 <9 2 u^ ! <9 2 w 0 2 ^ 2 T'| (At) 2 | ^5 3 t) d 2 w^ (ax) 2 

Lfe 2 VTt 2__a + + 

(1 + A 2 ) d A u 

12 d x* 


(ax) 2 > + 0[(At) 3 ] + 0[(At) 2 AX + 0[At(Ax) 2 ] + 0 [(ax) 3 


Assuming Eqs. (2.51) and (2.55), (e 2 )™ and (e 3 )™ reduce to 
( d 3 u\ n [(ax) 2 + a 2 (At) 2 ] 


(j,n)G 0 (2.68) 


\dx 3 ) . 


W = -(«!» 

and 


12 


+ 0[(At) 3 ] +0[(At) 2 Ax] +0 [a1(ax) 2 ] + 0[(ax) 3 ] (j,ri)£Cl (2.69) 


(e 3 )" = - ] + e>[(At) 3 ]+0[(At) 2 Ax]+0[Af(Ax) 2 ]+C>[(Ax) 3 ] (j, n) G 0 (2.70) 

respectively. 

Hereafter, for any u, let the system of equations defined by Eqs. (2.14), (2.15), and (2.59) be referred to 
as the basic form of the a(3) scheme while that defined by Eqs. (2.61)-(2.63) be referred to as the forward 
marching form of the a(3) scheme. Because (i) Eqs. (2.14), (2.15), and (2.37) <G> Eqs. (2.42)-(2.44) if A ^ 0, 
and (ii) Eqs. (2.59)-(2.63) are special cases of Eqs. (2.37), and (2.41)-(2.44), respectively, the basic form of 
the a(3) scheme its forward marching form. Thus the essential conditions represented by these or other 
equivalent forms may be referred to simply as the a (3) scheme. 

With the above definitions, the expressions on the right sides of Eqs. (2.64), (2.69) and (2.70) represent 
the truncation errors of Eqs. (2.61)-(2.63), respectively, if the forward marching form of the a(3) scheme is 
considered as a solver of the system of PDEs Eqs. (1.1) and (2.49). According to Eqs. (2.69) and (2.70), 
(e 2 )” — » 0 and (e 3 )" — > 0 as At, ax —> 0, regardless how At and ax are related when At, Ax — > 0. On 
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the other hand, Eq. (2.64) implies that (ei)" 
subjected to the condition 

(ax ) 4 

At, 


0 as At, ax — > 0 only if the mesh refinement procedure is 
0 as At, Ax — > 0 (2-71) 


Thus the a( 3) scheme is consistent with the system of PDEs Eqs. (1.1) and (2.49) if and only if Eq. (2.71) 
is satisfied. 

At this juncture, we offer the following remarks: 

(a) Let At/ ax be held as constant as At, Ax — > 0. Then for this mesh refinement procedure, Eqs. (2.64), 
(2.69), and (2.70) imply that the truncation errors for Eqs. (2.61)-(2.63), respectively, are third order, 
second order, and second order in At, and ax. 

(b) Because (i) each of the two decoupled subsystems in each of Eqs. (2.14) and (2.15) is STI invariant by 
itself, and (ii) the system Eq. (2.37) is also STI invariant if a and (3 are parameters independent of 
( j,n ), by the definition of STI invariance one can easily see that the basic form of the a(3) scheme is 
STI invariant. 

(c) Let q(j, n) = q 0 (j, n), ( j , n) € fl, be a solution to the basic form. Then, by substituting q[j, n) = q 0 (j , n) 
into the basic form, one obtains a system of identities involving q 0 (j,n), ( j,n ) € Q. Due to the STI 
invariance of the basic form, the above system of identities is equivalent to that obtained by substituting 
q(j,n) = Uq 0 (—j,—n) into the basic form. As such q(j,n) = q 0 (j,n), (/, n) G ft, represent a solution 
to the basic form •<=>■ q(j,n) = Uq 0 (~j,—n), ( j,n ) € O, represent another solution to the basic form. In 
other words, the STI image of a solution to the basic form is also a solution and vice versa. Obviously 
this conclusion is valid for other STI invariant forms of the a(3) scheme. 

Next, the forward marching form Eqs. (2.61)-(2.62) will be cast into a matrix form. Let 


\ def 

Co(^) = 


1 

—v 

(1 + 2i/ 2 )/3 


- ( \ def 

, c+(u) = 


1 

-(1 + A) | , S_(u) 

(2/3) (1 + v + v 2 ) 


def 


l 

i-i/ 

(2/3) (1 — v + v 2 ) 


(2.72) 


-(1 + v)/2 \ 




def 


do{v) = \2v , d+{v) (1 — 2i/)/2 , d-{v) = -(l + 2i/)/2 

-3/ V (3/ 2 ) / V (3/2) ) 

2 -2v (2/3) (1 + 2v 2 ) 

Qo{v) "= d 0 (v) [c 0 (iy)Y = | 2v —2v 2 (2/3)i/(l + 2v 2 ) 

-3 3i/ — (1 + 2v 2 ) 

' —(1 + v)/2 (1 + v) 2 /2 — (1 + i/)(l + v + v 2 )/?>' 

Q+{v) d+(v) [c + (is)} 1 = I (1 — 2v)/2 —(l — 2v)(l + v)/2 (1 — 2v)(l + v + i/ 2 )/3 


def 


3/2 


(3/2) (1 + v) 


l + v + v 1 


and 


(2.73) 


(2.74) 


(2.75) 


— (1 — v)/2 -(X-vf/2 -{l-v){l-v + v 2 )/Z 

Q.( y v) =d-{v)[c-{v)] t =\-{l + 2v)/2 — (1 + 2i/)(l — v)/2 -(1 + 2i/)(l - v + ^ 2 )/3 (2.76) 


3/2 


(3/2) (1 — v) 


1 — v + v 2 


Hereafter c* denote the transpose of any column or row matrix c. By using Eqs. (2.31) and (2.74)-(2.76), 
the forward marching form can be cast into the matrix form: 


q(j,n) = Qo{v)q{j,n- 1) + Q + {v)q(j + l,n- 1) + Q-{v)q{j -l ,n- 1), (j,n) £ fl (2.77) 

Here the reader is warned that the notations <5+(;/) and Q-(u) used in earlier CESE papers are now replaced 
by Q~(v) and Q + (v), respectively. As such, the terms Q-(u)q(j+ 1, n—T) and Q + {v)(j— 1, n—T) in Eq. (3.48) 
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of [71] appear here as Q+(v)q(j + 1 ,n— 1) and Q_(v)(j — 1 ,n — 1), respectively. Also note that each of 
Q o(^), Q+(v), and Q_(v) is in the form of dc 1 where c and d are 3x1 column matrix. Thus each is a matrix 
of rank one (see pp. 80-82 in [74]). Rank-one matrices are singular and have many interesting properties. 
As an example, the eigenvalues of Qo(v) are 0, 0, and \co(v)f do (is) with do ( is) being the eigenvector of the 
last eigenvalue. 

To facilitate the proof of the STI invariance of the forward marching form, first we will introduce some 
basic concept. Note that, for any set of variables xg, yg, t = 1, 2, the conditions 

xi + yi = x 2 - V 2 and x\ - yi = x 2 + 2/2 (2.78) 

Xi = x 2 and y\ = — j/ 2 (2.79) 

Thus, the image of Eq. (2.78) under any one-to-one mapping 

(xg,yg) <-> (x'gyy'g), f = l,2 (2.80) 

i.e., 

x'i + y[ = x' 2 - y' 2 and x\ - y[ = x' 2 + y 2 (2.81) 

<t=> the image of Eq. (2.79) under the same mapping, i.e., 

x'i = x 2 anc l Vi — ~V 2 (2.82) 

where the variables x\ and y], £ = 1,2, may or may not be related to xg,yg, l = 1,2. Moreover, in case 
that these two sets of variables are related, the condition Eq. (2.78) (or its equivalent Eq. (2.79)) may or 
may not be equivalent to the condition Eq. (2.81) (or its equivalent Eq. (2.82)). If the mapping Eq. (2.80) 
is such that Eq. (2.78) <t=> the image under this mapping (i.e., Eq. (2.81)), then Eq. (2.79) (the equivalent of 
Eq. (2.78)) <t=> Eq. (2.82) (the equivalent of Eq. (2.81)). Eq. (2.80) with x\ = xg and y' e = yg, l = 1,2, is an 
example of such mapping while Eq. (2.80) with x = ye and y] = xg, £=1,2, is not. 

To prove the STI invariance of the forward marching form, Note that: (i) the basic form of the a(3) 

scheme its forward marching form for any choice of q(j , n), ( j , n) € f2; and (ii) the STI images of the basic 

and forward marching forms, respectively, are obtained from the basic and forward marching forms through 
the mapping Eq. (2.30), i.e., through replacing q(j, n) in the basic form and the forward marching form with 
Uq(—j,—n), ( j,n ) £ O. From the above observations and the illustration given in the last paragraph, one 
concludes that the STI image of the basic form <£=> that of the forward marching form. Because the basic 
form is STI invariant, i.e., the STI image of the basic form <=> the basic form itself, Now we arrive at the 
conclusion that the forward marching form the basic form <=>■ the STI image of the basic form the STI 
image of the forward marching form. Thus the forward marching form its STI image, i.e., the forward 
marching form is STI invariant. QED. 

With the above preliminaries, the backward marching form of the a(3) scheme will be developed in 
Sec. 2.6. 

2.6. The backward marching forms of the a(3) scheme 

The STI invariance of the forward marching form of the a(3) scheme implies that Eq. (2.77) <t=> its STI 
image, i.e., 

Uq(—j, —n) = Q 0 (v)Uq(-j,-n+l) + Q + (v)Uq(-j - 1, -n+ 1) + Q-(is)Uq(-j + l,-n + 1), (j,n) G O 

(2.83) 

Moreover, by multiplying Eq. (2.83) from left using the matrix U and using Eq. (2.33), one concludes that 
Eq. (2.83) 

q(-J,-n) = Qo(v)q(~j, -n+l) + Q-(is)q(-j-l,-n + l) + Q+(v)q(-j + l,-n + l), ( j,n ) GO (2.84) 
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where 


QoM = f UQ 0 (u)U = 


( — (1 + z/ )/2 

Q-M d = UQ+{v)U = — (1 — 2is)/2 

V 3/2 

and 


/ 2 2i/ (2/3) (1 + 2zA) \ 

-2i/ — 2z/ 1 2 * -(2/3+(l + 2z+ 

\ -3 -3i/ -(l + 2z+ / 


-(l + ^) 2 /2 
-(1-2i/)(1 + i/)/ 2 
(3/2) (1 + i/) 


— (1 + i/)(l + v + z+/3 

— (1 — 2z/)(l + i/ + zA) /3 

1 + i/ + z/ 2 


/-(l-i/)/2 
Q+(v) = UQ_(is)U=\ (1 + 2z/)/2 
V 3/2 


(1-z,) 2 /2 
-(1 + 2i/)(1-z/)/2 
(3/2) (1 — z/) 


— (1 — i/)(l — z' + zA)/3 
(1 + 2z/)(l — z/ + z/ 2 )/3 
1 — z/ + z/ 2 


(2.85) 


(2.86) 


(2.87) 


By replacing the “dummy” indices —j and —n everywhere in Eq. (2.84) with j and n, respectively, one can 
see that the system Eq. (2.84) is identical to the system 

q{j, n) = Qo(v)q{j, n + 1) + Q+{v)q(j + 1, n + 1) + Q-(u)q(j - 1, n + 1), (j, n) G Cl (2.88) 

Because the mesh variables at O’, n) can be determined in terms of those at (j — 1 , n + 1), ( j,n + 1), and 
(j + 1 , n + 1) using Eq. (2.88), hereafter Eq. (2.88) (which is equivalent to other forms of the a( 3) scheme) 
will be referred to as the backward marching form of the a(3) scheme. 

According to Eqs. (2.74) and (2.85), Qo(z') = Udo{v) [co(i')] t U. Because Udo{v) and [c'o(z')] 4 U are 3x1 
column matrix and 1x3 row matrix, respectively, Q o(z') is a rank-one matrix. Similarly, Q_(z/) and Q+(z^) 
are also rank-one matrices. 

Eq. (2.88) was derived using the STI invariance of the forward marching form of the a(3) scheme. 
Alternatively, it can also be derived from the basic form. To proceed, note that: (i) by replacing the indices 
j and n everywhere in Eq. (2.14) with j + 1 and n + 1 and using the fact that (j, n) G Cl <i+ (j — 1, n — 1) G 
^ ^ 0 + A n + 1) £ ^ one can see that the system Eq. (2.14) is identical to the system 


, \ 2(1 — v + v 2 ) 

U + (1 - v)Ux H Uxx 


„ , 2(1 — z/ + z/ 2 ) in+1 

u - (1 - v)Ux H Uxx 


(j,n)€Cl (2.89) 


3 + 1 


(ii) by replacing the indices j and n everywhere in Eq. (2.15) with j — 1 and n + 1 and using the fact that 
(j, n) € Cl <+> () + l, n — 1) G f 1 <+> (j — 1, n + 1) € Cl, one can see that the system Eq. (2.15) is identical to 
the system 


n , >1 , 2(1 + is + v 2 ) 

U - (1 + V)Ux H 5 Uxx 


, . 2(1 + v + v 2 ) 

u + (1 + v)Ux H 5 u S x 


n+1 


(),n)GO (2.90) 


1 


and (iii) by replacing the index n everywhere in Eq. (2.59) with n+1 and using the fact that (j. n) G Cl •+> 

(j, n — 1) G Cl <+> (j, n + 1) G Cl, one can see that the system Eq. (2.59) is identical to the system 


1 + 2zA 

U - VUx H Uxx 


J 3 


1 + 2z/ 2 

« + vux H Wji 


n+1 


0, n) G fi 


(2.91) 


J j 


As such the system Eqs. (2.89)~(2.91) are identical to Eqs. (2.14), (2.15), (2.59), respectively. 

For each ( j,n ) G Cl, Eqs. (2.14), (2.15), and (2.59) form a linear system of three equations for the three 
mesh variables n", ( u s )", and (n®®)". Eqs. (2.89)-(2.91) form another system. Moreover, one can see that, 
under the mesh variable mapping 


q(j, n) <-> C7 j, n) , g(j, n - 1) <-> g(j, n + 1) , 

q{j + l,n - 1) <->■ Uq(j - l,n + 1), and gO - 1, n - 1) <-> g(j + 1, n + 1) 
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Eqs. (2.89)-(2.91), respectively, are the images of Eqs. (2.14), (2.15), and (2.59) and vice versa. By using 
the concept introduced earlier in a discussion involving Eqs. (2.78)-(2.82), one concludes that the solution 
to Eqs. (2.89)-(2.91) must be the image of Eq. (2.77) (i.e., the solution to Eqs. (2.14), (2.15) and (2.59)) 
under the same mapping. In other words, the solution to Eqs. (2.89)-(2.91) is 

Uq{j,n ) = Q 0 (is)Uq(j,n+ 1) + Q + {v)Uq{j - l,n + 1) + Q-{v)Uq{j + l,n + 1), (j,n) G (2.93) 

By multiplying Eq. (2.93) from left using the matrix U and using Eqs. (2.33) and (2.85)-(2.87), one has 
Eq. (2.88). QED. 

As a preliminary for the developments in Sec. 3, in the following, important algebraic relations involving 
Qo{v), Q+(v), <3_(A), Q o(i'), Q+(v), and Q-{v) will be extracted from the STI invariance of the a(3) scheme. 

2.7. Algebraic relations associated with STI invariance 

Let (jo, n 0 ) G 12 be any given fixed mesh point. Let q(j 0 , n 0 ), q(j 0 ± 1, n G ), and q(j 0 ± 2, n 0 ), respectively, 
be the arbitrary initial data specified at (j D , n 0 ), (j G ± 1, n 0 ), and (j Q ± 2, n D ), respectively. Let q(j 0 , n 0 + 1), 
and q(j 0 ± l,n Q + 1) be specified in terms of the mesh variables at the n Q th time level using the forward 


marching form Eq. (2.77), i.e., 

q{jo, n a + 1) = Qo(v) q{jo, n a ) + Q+(v) q{j Q + 1, n Q ) + Q-(v) q{jo ~ 1, n a ) (2.94) 

q(j 0 + 1, n a + 1) = Qo(^) q{jo + 1, n Q ) + Q+{v) q{j a + 2, n Q ) + Q-(^) q{j a , n 0 ) (2.95) 

and 

q(jo - 1, n a + 1) = Q 0 {v) q{j Q - 1, n Q ) + <5+(^) q(j a , n 0 ) + Q-{v) q(j a - 2, n Q ) (2.96) 

On the other hand, because Eq. (2.77) <t=> Eq. (2.88), q{j 0 ,n 0 + 1), q(j a ± l,n 0 + 1), and q(j 0 ,n 0 ) must also 

be linked by Eq. (2.88), i.e., 

q(jo, n-o) = Qo(^) q{jo, n Q + 1) + Q + {v) q{j a + 1, n 0 + 1) + Q-{v) q{j a -l,n 0 + 1) (2.97) 

Substituting Eqs. (2.94)-(2.96) into (2.97), one has 

[Qo(v)Qo(v) + Q+(v)Q-(v) + Q-{v)Q+{v) - l\ q{jo, n 0 ) 

+ [Qo{v)Q + (v) + Q+(v)Qo(v)]q(jo + l,n 0 ) + [Qo{v)Q_{v) + Q-{v)Q 0 {v)\q{jo - l,n 0 ) (2.98) 

+ Q+(y)Q+(y) q(j 0 + 2, n Q ) + Q-(y)Q_(y) q(j a - 2 ,n D ) = 0 


where I is the 3x3 identity matrix and 0 is the 3x1 null column matrix. 

Because Eq. (2.98) must be valid for any choice of q(j a ,n 0 ), q(j 0 ± l,n 0 ), and q(j 0 ± 2 ,n 0 ), the coeffi- 
cients matrices in front of these column matrices must be null identically, i.e., 


Qoiy)Q 0 (y) + Q + (v)Q_(y) + Q-(y)Q + (y ) = / (2.99) 

Q 0 {u)Q+{v) + Q+{v)Q 0 {v) = 0 (2.100) 

Q 0 {v)Q-{v) + Q-{y)Q 0 {v) = 0 (2.101) 

Q+(v)Q+(v) = 0 (2.102) 

and 

Q-{y)Q-{v) = 0 (2.103) 


where 0 is the 3x3 null matrix. As an example, one can prove Eq. (2.99) by substituting into Eq. (2.98) 
each of the following sets of the initial data: (i) q(j Q ± 1 ,n Q ) = q{jo ± 2 ,n D ) = 0 and q{j 0 ,n 0 ) = (1,0,0)*, 
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(ii) q{j 0 ± 1, n 0 ) = q(j a ± 2, n a ) = 0 and q(j a , n Q ) = (0, 1, 0)*, and (iii) q{j a ± 1, n 0 ) = q[j a ± 2, n D ) = 0 and 

q{jo,no) = ( 0 , 0 , 1 )*. 

Similarly, by substituting the backward marching relations 

n a - 1) = Qo(v) q{jo, n a ) + Q+(v) q{j Q + 1, n Q ) + Q-{v) q{j Q - 1, n 0 ) (2.104) 

q(jo + 1, n Q - 1) = Qo(^) q{jo + 1, n 0 ) + Q+(^) q{jo + 2, n Q ) + Q-{v) q(j a , n 0 ) (2.105) 

and 

q(jo - 1, n Q - 1) = QoM - 1, n 0 ) + Q+(v) q(j 0 , n 0 ) + Q-{v) q(j Q - 2, n Q ) (2.106) 

into the forward marching relation 

q{jo, n a ) = Q 0 {v) q(j 0 , n Q - 1) + Q+{v) q{j Q + 1, n c - 1) + Q-(y) q{j Q ~ 1, n a - 1) (2.107) 

one has 

[Qo{v)Qo{v) + Q+(v)Q-{v) + Q_{v)Q + (v) - l]q(j 0 ,n 0 ) 

+ [Qo(v)Q+(v) + Q + (v)Q 0 {v)]q{j o + l,n G ) + [Q 0 (is)Q-(v) + Q-(v)Q 0 (is)]q(j o - l,n 0 ) (2.108) 

+ Q+{v)Q+(v) q{jo + 2, n 0 ) + Q-(y)Q-(y) q(j Q - 2 ,n Q ) = 0 

Because Eq. (2.107) must be valid for any choice of q(jo,no), q(jo i l,n G ), and q(j Q ± 2 ,n 0 ), one concludes 
that 

QoMQoM + Q+(v)Q-(v) + Q-(y)Q+(y) = / (2.109) 

Qo(y)Q + {v) + Q + {y)Q 0 {y) = 0 (2.110) 

<5o(^)<3-(^) + <5-(^)Qo(^) = 0 (2.111) 

Q+(i/)0+(i/) = O (2.112) 

and 

Q^(y)Q-{v) = 0 (2.113) 

By using Eqs. (2.32) and (2.85)-(2.87), it can be shown that: (i) Eq. (2.99) <t=> Eq. (2.109) <t=> 

Qo{v)UQq{v) + Q-{v)UQ-{v) + Q + {v)UQ + {v) = U (2.114) 

(ii) Eq. (2.100) Eq. (2.111) 

Qo{v)UQ+{v) + Q-(v)UQo{v) = 0 (2.115) 

(iii) Eq. (2.101) Eq. (2.110) 

Q 0 {y)UQ-[y) + Q+[y)UQ 0 {y) = 0 (2.116) 

(iv) Eq. (2.102) Eq. (2.113) 

Q-{y)UQ+{v)=Q (2.117) 

and (v) Eq. (2.103) Eq. (2.112) 

Q + (y)UQ-{y) = 0 (2.118) 

2.8. Other invariant properties and related algebraic relations 

By using Eqs. (2.32) and (2.74)-(2.76), one can show that 

Qo{-v) = UQ 0 (v)U, Q-{-v) = UQ+(is)U, and Q + (-u) = UQ-{v)U (2.119) 
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By using Eqs. (2.85)-(2.87), one can also show that Eq. (2.119) <t=> 

Q 0 {-v) = UQ 0 {v)U, Q-(-v) = UQ+(v)U, and Q+(-i/) = UQ-{v)U (2.120) 

As will be shown, the above relations are linked with other invariant properties of the a(3) scheme. 

Let the advection speed a in Eq. (1.1) be considered as a variable parameter. Let u = u(x,t\a) be a 


solution to Eq. (1.1), in the domain — oo < x,t,a < + 00 , i.e., 



du(x , t; a) du(x, t\ a) 

dt +a dx 

— 00 < x,t,a < +00 

(2.121) 

Let 

x' = f —x, t' = f t, and a' = —a, 

—00 < x,t,a < +00 

(2.122) 

and 

u(x,t;a ) = u(—x,t;~ 

-a) 

(2.123) 

Then (i) Eq. (2.121) 



du(x' , t'\ a') ,du(x', t’\ a!) n 

dt' +a dx' 

—00 < x ’ , t', a' < +00 

(2.124) 

and (ii) 

d d d 

di’ = dt an d? = 

d 

dx 

(2.125) 

Thus one concludes that Eq. (2.121) <t=> 



du(x,t;a) du(x,t;a) n 

dt +a dx ~ ’ 

—00 < x, t < +00 

(2.126) 


In other words, if u = u{x,t\a ) is a solution to Eq. (1.1), so must be u = u(x,t\a) and vice versa. Because 
the one-to-one mapping 


(x, t, a) <— > (— x, t, — a), —oo<x,t,a<+oo (2.127) 

represents a combined spatial-reflection (parity) and advection direction reversal (ADR) operation, hereafter 
(i) a pair of functions such as u and u will be referred to as the PADR images of each other; and (ii) a PDE 
such as Eq. (1.1) is said to be PADR invariant if the PADR image of a solution is also a solution and vice 
versa. 

Because v = aAt/Ax, the numerical analogue of Eq. (2.127) is 

(j, n) (— j,n ) and v <-> — v (2.128) 

Motivated by an argument similar to that leads to Eq. (2.30) for STI mapping, the PADR mapping for the 
a (3) scheme is defined by 


q(j,ri) Uq(—j,ri) and v — u, (2.129) 

Thus the PADR image Eq. (2.77) is 

Uq[—j, n) = Qo(—v)Uq(—j,n — 1) + Q + {-v)Uq(-j - l,n - 1) + Q-(-v)Uq(-j + 1 ,n- 1), (j,n) G Q 

(2.130) 
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By using Eqs. (2.32) and (2.119), it can be shown that Eq. (2.130) <t=> 

n ) = Qo(v)q(-j,n - l) + Q-(v)q(-j - l,n - l) + Q+(v)q(-j + l,n - l), ( j,n ) G fl (2.131) 

By replacing the dummy index —j with j everywhere in Eq. (2.131) and using the fact that (— j,n ) G 12 
(j,n) G fi, one concludes that Eq. (2.131) <t=> Eq. (2.77). Thus Eq. (2.77) is PADR invariant, i.e., it is 
equivalent to its PADR image. 

By exchanging the roles of x and t, one can define invariance under a combined time reversal and 
advection direction reversal operation. Because (i) this operation is equivalent to a STI operation followed 
by a PADR operation or vice versa, and (ii) Eq. (1.1) and the a(3) scheme are invariant under both STI 
and PADR operations, one concludes that Eq. (1.1) and the a(3) scheme are also invariant under the new 
operation. In fact, invariance of the a(3) scheme under this new operation can be proved using Eq. (2.120) 
(which is equivalent to Eq. (2.119)). 

3. von Neumann analysis 

Let G{y, 9) be a 3 x 3 nonsingular complex matrix function of v and the phase angle 9 such that 

q(j, n ) = e y ' e [G(v, 0)]" 6, ( j , n) G fi; — oo < v, 9 < +oo ( i = \/ — 1) (3.1) 

is a solution to Eq. (2.77) for all possible complex constant 3x1 column matrices b. (Note: because 

[G(v, 9)] n = f |[G(i/, 0)] -1 j for an integer n < 0, [G(v, 9)] n is not defined if n < 0 unless [G(v, 0)] _1 exists, 
i.e., G(v,9) is nonsingular.) By substituting Eq. (3.1) into Eq. (2.77), one has 

[G(i/, 9) - Qo(v) - e l0 Q+{v) - e~ ie Q-{v)] [G{y, 9)} n 6 = 0, n = 0, ±1, ±2, . . . (3.2) 

Because (i) [G(y, 0)]° = /, and (ii) 6 can be any complex constant 3x1 column matrix, Eq. (3.2) <£=> 

G(v, 9) = Q 0 {y) + e l6 Q+(u) + e~ i0 Q.{u) (3.3) 

By definition, G(y,9) is the amplification matrix of the forward marching form of the a(3) scheme. Because 
Qo{v), Q+(v), and Q~{v) are real matrices, Eq. (3.3) implies that 


G{is,-9) = G(v,9) (3.4) 

Hereafter M denotes the complex conjugate of any matrix M. Also, with the aid of Eq. (2.119) and the 
relation U = G -1 , one has 

G(—u, 9) = UG{v, ~9)U = UG{v, -9)U~ l (3.5) 

In the following, we will show that the a(3) scheme must be neutrally stable when it is stable. 

3.1. Neutral stability of the a( 3) scheme 

By using Eqs. (2.114)-(2.118), one can show easily that 


U[Q 0 (v) + e i0 Q-{v) + e- l0 Q+{v)\ U [Q 0 [y) + e l0 Q+{v) + e~ i0 Q-(is)] 

= [Qo(v) + e i0 Q+(v ) + e~ i0 Q-(v)] U [Q 0 (v) + e l0 Q.(u) + e~ i0 Q+(u)] U = I 
Thus G(y,9) defined in Eq. (3.3) is nonsingular and its inverse is 

G{y.9)\ 1 = U [Qo(v) + e l0 Q.(v) + e~ l0 Q + (y)} U 


(3.6) 


(3.7) 
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Indeed, with the aid of Eq. (2.85)-(2.87), Eq. (3.7) is what one obtains after substituting Eq. (3.1) into the 
backward marching form Eq. (2.88). Moreover, by using Eqs. (2.32), (3.3), (3.4), and (3.7), one has 

[G(i/, 0)] _1 = UG(v, 9)U _1 (3.8) 

For each (v,9), the three eigenvalues G(i/,6) will be denoted as <jz(v,9), £ = 1,2,3, and referred to as 
the amplification factors of the a(3) scheme. Because G(v, 9) is nonsingular, 

<7i(y, 9)^0, £=1,2,3; — oo < u, 9 < +oo (3.9) 

(see part (i) of Theorem 1 given below). Also, as will be shown, at(v,9), £ = 1, 2, 3, satisfy the following set 
condition: 


1 , 1 , 1 \ 

<7\(v, 9) ’ (T2(v, 9)’ cr 3 (iy,9)J 


{(?i(v,9), cr 2 (i/, 0), o- 3 (i/,6>)| , 


—00 < is, 9 < +00 


(3.10) 


Hereafter z denotes the complex conjugate of any complex number z. 

As a preliminary, first we introduce the following matrix theorems: 

Theorem 1 . Let A be a nonsingular N x N matrix with the eigenvalues A^, £ = 1 , 2 , . . . , N . Then (i) 
A# ^ 0, £ = 1, 2, . . . , TV; and (ii) the eigenvalues of A^ 1 are 1/Ae, £ = 1,2 ,,N. 

Theorem 2. Let A be a IV x N matrix with the eigenvalues At, £ = 1,2,..., iV. Then the eigenvalues 
of A, the complex conjugate of A, are At, £ = 1,2,..., TV. 

Theorem 3. Let A and B be two similar N x N matrices, i.e., there exists a nonsingular TV x N matrix 
S so that B = S' -1 AS 1 . Then A and B have the same eigenvalues, counting multiplicity. 

The proof of Theorems 1 and 2 is given in Appendix A while that of Theorem 3 is given on p. 45 of [76]. 

To prove Eq. (3.10), note that part (ii) of Theorem 1 implies that, for any (v, 9), the eigenvalues of 
[G(u, 0)] _1 are 1 /<jt(v,9), £ = 1,2,3. Next, by using Theorems 2 and 3, and the fact that (Z7 -1 ) -1 = U, 
one can see that the eigenvalues of the matrix on the right side of Eq. (3.8) are <Jt(y,9), £ = 1,2,3. Thus 
Eq. (3.10) now is an immediate result of Eq. (3.8). QED. 

An immediate result of Eq. (3.10) is 


111 
ai{v,9) <72 (v, 9) <7 3 {v,9) 


<Ji{v, 9) • cr 2 (v,9) • <r 3 (i/, 9) 


i.e., 


Wi{v,9)\ ' MM)] • W3(u,9)\ = 1 

For any given v, stability of the a(3) scheme requires that 

\<n{v,9)\<\, £=1,2,3 

Thus Eq. (3.11) implies that, for any given v, the a(3) scheme must be neutrally stable, i.e., 


Wt(y,9) | = 1, £=1,2,3 


(3.11) 


(3.12) 


(3.13) 


if it is stable. As such, Eq. (3.8) does not imply neutral stability of the a( 3) scheme. However, it does imply 
that the scheme can only be neutrally stable (i.e., non-dissipative) if it is stable. Here we have reached this 
conclusion without using the explicit form of at(v, 9), £=1,2, 3. 

At this juncture, note that one can obtain 


at(—v,9) = at(v,—9) = at(v,9), £=1,2,3 (3.14) 
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by using Eqs. (3.4) and (3.5) along with Theorems 2 and 3. 

Eq. (3.10) and (3.14) are the fundamental relations governing the eigenvalues of G(y , 0). In the following, 
we explore other properties of these eigenvalues. 

3.2. Characteristic equation of G(y,9) 

By using Eqs. (2.74)-(2.76) and (3.3), one has 


G[y,9) = 


^ 2 — cos 0 — ir/sin 0 

2i/(cos 0 — 1) + «(1 + v 2 ) sin 0 

2(1 + 2is 2 ) m 2iis(2 + z/ 2 ) . 

( 1 — cos 0) sin 0 

3 v ' 3 

2i/(l — cos 0) + i sin 0 

(2v 2 — 1) cos 0 — 2v 2 + iv sin 0 

2^(1 + 2v 2 ) 2i(l — v 2 ) . n 

v ’ (1 cos 0) + v ’ sin 0 

O O 

\ 3(cos0 — 1) 

3^(1 — cos 0) — 3 i sin 0 

2(1 + v 2 ) cos 0 — 1 — 2v 2 + 2 iv sin 0 

— 00 < 13,9 < +OO 




(3.15) 

It follows from Eq. (3.15) that (i) 

det [G{v,9)\ = — 1, —oo < v, 9 < +oo (3.16) 

and (ii) any eigenvalue cr of G(v, 9) must be a root of the characteristic equation: 

det [a I — G(v, 0)] = a 3 + h(v, 9)a 2 + h(v, 9)a +1 = 0 (3-17) 

where 

h(v, 9) =* —1 + dv 2 {\ — cos 9) — 2ivs\\\9 (3.18) 

The reader may be surprised by the simple result Eq. (3.16). However, by using Eq. (3.3) and the fact 
that each of Qo(v), Q+(v), and Q~(v) has the form dc * with c and d being 3x1 column vectors, an 
application of the fundamental definition of determinant (in which the Levi-Civita antisymmetric symbol is 
used) leads to the conclusion that det[G(+0)] must be independent of 0, i.e., det[G(+0)] = det[G(+0)]. As 
such, Eq. (3.16) now follows from the fact that G(v, 0) = U (see Eqs. (3.15) and (2.32)) and det(U) = — 1. 
Hereafter, for simplicity, the arguments v and 9 may be omitted if no confusion would arise. 

Because a i, 0 + and <73 are the eigenvalues of G, Eq. (3.17) implies that 

cr 3 + ha 2 + ha + 1 = (cr — cri)(cr — < 72 )(cr — 03 ) (3.19) 


for any complex variable cr. On the other hand, because Eq. (3.10) <+> 
{ai(v,0), a 2 (v,9), a 3 (v,6)} = ^ 


£Tl(+ 0) cr 2 (+0) cr 3 (l/, 0) 

1/oT, 1/ng, and l/o+ must also be the eigenvalues. Thus 


1 


a + ha + ha + 1 = [a — — cr — — cr — — 


cr 1 


1 


o' 2 


— 00 < v, 9 < +00 


(3.20) 


1 


<73 


(3.21) 


for any complex variable a. In the following, independently, it will be shown that indeed Eq. (3.19) implies 
Eq. (3.21). 
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Proof. Let a = 0. Then Eq. (3.19) implies that 


f 7 ! CT 2 <J 3 = —1 


(3.22) 


Eq. (3.9) is an immediate result of Eq. (3.22). Also, one can see that Eq. (3.21) -o- Eq. (3.22) if o = 0. 
Let <7^0. Then, by replacing a with l/a in Eq. (3.19), one has 


1 h 


=3+=2+ = + 1 — I A ~ a l ) I A “ °2 ) ( II 'As 


(3.23) 


Also, by using Eq. (3.22), one has 


cr 3 = (— ct/< 7 i)(— cr/< 7 2 )(— 0-/0-3) 


(3.24) 


Because the product of the expressions on the left sides of Eq. (3.23) and (3.24) equals to that on the right 
sides, we have 


cr 3 + her 2 + her + 1 = [a A] I a A) I a A 


<71 


1 


02 


1 


03 


(3.25) 


Eq. (3.21) is the complex conjugate form of Eq. (3.25). QED. 
Moreover, according to Eq. (3.18), 


h{—v , 9) = h(v, —9) = h(v , 9) 


(3.26) 


Thus Eq. (3.14) can also be derived directly from Eq. (3.17). 
In Sec. 3.3, we will prove the following proposition: 


\<re(y, 0)| = 1, £=1,2,3; — 00 < 9 < +00 if \u\ < 1/2 


(3.27) 


3.3. A proof of proposition Eq. (3.27) 

First we introduce the following well-established algebraic theorem: 

Theorem 4. Let ui, cr 2 , . . . , crjv' be the distinct roots of an algebraic equation of TVth order, i.e., 

cr N + aicr^ -1 + a 2 o N - 2 + . . . + ajv-icr + a,/v = 0 (3.28) 

where a 1 , a 2 , . . . , ajv are complex constant coefficients and cr is a complex variable. For each £=1,2,..., N', 
let me > 1 denote the multiplicity of the root cr^. Then 


N' 

y^me = N 

e=i 


(3.29) 


According to the above theorem, for any given (ir,9), the roots of the cubic equation Eq. (3.17) must 
fall into one of the following three mutually exclusive cases: (a) there is one triple root (multiplicity = 3); 
(b) there are one double root (multiplicity = 2) and one simple root (multiplicity = 1) and (c) there are 
three simple roots. 

Consider case (a). Then <j\ = cr 2 = cr 3 . Let cr Q denote the common value of o \ , cr 2 , and cr 3 . Then 
Eqs. (3.9) and (3.20) imply that (i) a 0 ^ 0 and (ii) 1/7TZ must also be a triple root of Eq. (3.17). Thus the 
only choice that will not contradict Theorem 4 is that cr 0 = 1 /oy, i.e., |er 0 | = |oi| = |cr 2 | = |cr 3 | = 1. 

Consider case (b). Without any loss of generality, one can assume a\ = cr 2 cr 3 . Again let a 0 denote 
the common value of and cr 2 . Then Eqs. (3.9) and (3.20) imply that (i) cr 0 0; (ii) cr 3 ^ 0; and (iii) 1/WZ 
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and 1/0-3 must also be a double root and a simple root of Eq. (3.17), respectively. Thus the only choice that 
will not contradict Theorem 4 is that a D = l/oy and 03 = 1/173, i.e., |er 0 | = |oi| = |£t 2 | = |cr 3 1 = 1. 

The conclusions reached above imply the following lemma: 

Lemma 1 . For any given ( v,6 ), the roots of Eq. (3.17) must all be of unit magnitude if any one of 
them is a multiple root. 

Thus, to prove Eq. (3.27), we need only to consider case (c), i.e., the case with 



ci yf cr 2 , er\ 7^ er 3 , 

and (72 yf <j 3 

(3.30) 

To proceed, each ere is expressed in 

its polar form, i.e., 




oe = ree x<t>e , 

£=1,2,3 

(3.31) 

where, because of Eq. (3.9) 

r e = \ae\ > 0, 

£=1,2,3 

(3.32) 

Moreover, we assume that 

-7r < fa < 7 r, 

£=1,2,3 

(3.33) 

Thus, for each ere, the value of e % ^ 1 

is unique. It follows from Eqs. (3.31) and (3.32) that 



1/ct = (1 /re)e t4 ‘ t , 

£=1,2,3 

(3.34) 


Also, by using Eqs. (3.31)-(3.33), Eqs. (3.30) can be expressed as the following ordered pair inequalities: 

(n,fa) 7^ (r 2 ,<fa), (ri,fa) ^ (r 3 ,fa), and (r 2 , fa) ± (r 3 > fa) (3-35) 

The distribution of fa, <fi 2 , and <j > 3 must fall into one of the following mutually exclusive cases: (cl) all 
have distinct values; (c2) two of them have the same value while the third assumes a different value; and 
(c3) all have the same values. In the following, these sub-cases will be discussed separately. 

Consider case (cl) where 

fa 7^ fa, fa 7^ fa, and (j) 2 7^ fa (3.36) 

Because Eqs. (3.20) and (3.34) imply that (1 /re)e l<t>t , £ = 1,2,3, must also be roots of Eq. (3.17), Eq. (3.36) 
implies that the only choice that will not contradict Theorem 4 is that re = 1/re, be., re = 1, £ = 1,2,3. 
Thus, for case (cl), again we have \ai\ = |cr 2 1 = |cr 3 | = 1. 

Consider case (c2) where, without any loss of generality, one can assume that 

fa = fa^ fa (3.37) 

Because of (3.35), Eq. (3.37) implies that 

n 7^ r 2 (3.38) 

By using Eqs. (3.37) and (3.38) along with the fact that (1 /re)e l<l>i , £=1,2, 3, must also be roots of Eq. (3.17), 
one concludes that the only choice that will not contradict Theorem 4 is that r\r 2 = 1, r\ ^ 1, r 2 ^ 1, and 
r 3 = 1. Thus, for case (c2), (i) one of the roots is of unit magnitude while the other two are not; and (ii) 
the product of the magnitudes of the two roots which are not of unit magnitude is one. 

Consider case (c3) where 

fa = fa = fa (3.39) 

Because of Eq. (3.35), Eq. (3.39) implies that 

r\ yf r 2 , r\ r 3 , and r 2 ^ r 3 (3.40) 
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By using an argument similar to that invoked in the discussion of case (c2), one concludes that, for case 
(c3), again (i) one of the roots is of unit magnitude while the other two are not; and (ii) the product of the 
magnitudes of the two roots which are not of unit magnitude is one. 

As a result of the above discussions, we have the following lemma: 

Lemma 2. For any given (v,9), the case with at least one of the roots of Eq. (3.17) not being of unit 
magnitude may occur only if it meets the following conditions: (i) one and only one of ri, V 2 , and r 3 is of 
unit magnitude; and (ii) the two roots that are not of unit magnitude share the same phase angle and the 
product of their magnitudes is one. 

Consider any case that meets the conditions referred to in Lemma 2. Then, without any loss of generality, 
one may assume that 

nr 2 = 1, r\± 1, r 3 = 1, and fa = fa = (3-41) 

where (f> denotes the common value of fa and fa. Moreover, by using Eqs. (3.31) and (3.34), Eqs. (3.18), 
(3.19), and (3.21) imply that 


nr 2 r 3 e i ^ 1+cl>2+ ^ — ^ c *( 0 1+02+03) _ _i 

rir 2 r 3 

nr 2 e i ^ 1+(/>2 ^ + r\r 3 e i ^ 1+ ^ + r 2 r 3 e i ^ <t>2+ ^ = ^ c *( 0 1+02) _(_ ^ c »Oi+<fe) _|_ ^ c *( 02 + 03 ) 

'V'j nr 3 r 2 r 3 

= — 1 + 4j/ 2 (1 — cos 6) + 2ii/ sin 0 

and 

+ r 2 e 1 ^ 2 + r 3 e l< ^ 3 = — e^ 1 H e*^ 2 H e l ^ 3 = 1 — 4*/ 2 (l — cos0) + 2**/sin0 

n r 2 r 3 


rie 


,*0i 


Because of Eq. (3.32), Eq. (3.42) 
and 


r 2 


rs 




nr 2 r 3 = 1 


gi(01+02+03) _ 

By using Eqs. (3.45) and (3.46), Eq.(3.43) 


(3.42) 

(3.43) 

(3.44) 

(3.45) 

(3.46) 


1 


ne *^’ 1 + r 2 e + r 3 e l<t>3 = — e t<t ‘ 1 

r 1 


<i r 2 i 3 

Because Eq. (3.47) is the complex conjugate of Eq. (3.44), we conclude that Eqs. (3.44)-(3.46) <t=> Eqs. (3.42) 
(3.44), i.e., Eqs. (3.44)-(3.46) represent all the independent constraints imposed on ry and <pe, l = 1,2,3. 
Note that Eq. (3.45) is satisfied by Eq. (3.41) automatically. However, Eqs. (3.41) and (3.46) imply tha 


e *^ 2 H e = 1 — 4 z/ 2 (1 — cos 9) — 2*j/sin 6 (3.47) 

r 3 


1 on re and fa, t — 1,2, 3. 
(3.41) and (3.46) imply that 

g*03 _ g — *(01+02) _ g — 2*0 


Let 


def 

p = r 1 


r — ' L 

Then, with the aid of Eqs. (3.41) and (3.48), Eq. (3.44) can be expressed 


as 


where 
Eq. (3.50) 


^ — e 2 *0 = 1 — 4 j/ 2 (1 — cos0) + 2*z/sin0, 

p > 0 and p ^ 1 


— 7r < </> < 7r; p > 0 and p ^ 1 


/M = - ■ 1 


P 4 — ) 
P 


(3.48) 

(3.49) 

(3.50) 

(3.51) 


/(p) cos </> — cos(2 fa = 1 — 4i* 2 (1 
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and 


f(p) sin <j> + sin(2</>) = 2^sin0, — n < (j> < 7r; p > 0 and 1 (3.53) 

Thus, given any ( is, 9 ), Eqs. (3.52) and (3.53) must admit a solution for p and 4> in the specified domain if 


the case Eq. (3.41) indeed exists. 

To explore Eqs. (3.52) and (3.53), note that (i) 

[. f(p ) cos<f) — cos(20)] 2 + [f(p) sin</> + sin(2^)] 2 = [l — 4 i^ 2 (1 — cos0)] 2 + [2^sin0] 2 (3.54) 

is a direct result of Eqs. (3.52) and (3.53); (ii) 

[f(p) cos (j) - cos(20)] 2 + [f(p) sin (j> + sin(20)] 2 = [f(p) - l] 2 + 2 f(p) [1 - cos(3</>)] (3.55) 

and (iii) 

[l — Av 2 (l — cos0)] 2 + [2issm9] 2 = 1 — 4is 2 (1 — 4is 2 )(l — cos 9) 2 (3.56) 

Next, because (i) the minimum of f{p) in the domain p > 0 occurs at p = 1 and (ii) /( 1) = 2, we have 

f(p) >2 if p > 0 and p ^ 1 (3.57) 

Combining Eqs. (3.55) and (3.57), and using the fact that 1 — cos(3 <j>) > 0 for all (f>, one has 

[ f(p ) cos </> — cos(2</>)] 2 + [f(p) sin <j> + sin(2(/>)] 2 >1 if p > 0 and p ^ 1 (3.58) 

On the other hand, because 1 — 4^ 2 > 0 if \is\ < 1/2, Eq. (3.56) implies that 

[l — 4v 2 (\ — cos$)]~ + [2j/sin0] 2 < 1 if \v\ < 1/2 (3.59) 

Combining Eqs. (3.58) and (3.59), one arrives at the conclusion that, for all 6 , 

[f(p) cos (j) — cos(2^)] 2 + [f(p) sin</> + sin(2^)] 2 > [l — 4 i^ 2 (1 — cos0)] 2 + [2^sin0] 2 (3.60) 


i.e., Eq. (3.54) cannot be satisfied, if (i) \v\ < 1/2; and (ii) p > 0 and p / 1. Because Eq. (3.54) is a direct 
result of Eqs. (3.52) and (3.53), this implies that, for any 9 , Eqs. (3.52) and (3.53) admit no solution for p 
and </> in the specified domain, i.e., the case Eq. (3.41) does not exist. In turn, this implies that, for all 9 , 
the roots of Eq. (3.17) are all of unit magnitude if \v\ < 1/2, i.e., Eq. (3.27) has been proved. QED. 

Note that, by itself, Eq. (3.27) does not imply that the a(3) scheme is stable when \u\ < 1/2. According 
to the discussion given in Sec. 4 of [1] and Secs. 3 and 4 of [64], the fact that all its amplification factors are 
of unit magnitude for all phase angles at a value of v insures that the a(3) scheme is stable at this particular 
value of v only if the amplification matrix G(v,9) is nondefective for all 9. By definition [76], for any given 
(is, 9), the 3x3 matrix G(is, 9) is nondefective if the dimension of its eigenspace (i.e., the vector space spanned 
by its eigenvectors) is three. On the other hand, G(is, 9) is defective if the dimension of its eigenspace is less 
than three. It can be shown that G(is, 9) is defective if and only if its Jordan canonical form [76] contains 
a Jordan block with its dimension > 1. As such G(v,9) is defective if and only if (i) it has an eigenvalue 
with multiplicity > 1 and (ii) the Jordan block associated with the eigenvalue has a dimension > 1. For the 
defective case, stability requires a stronger condition, i.e., the magnitude of the eigenvalue associated with 
the Jordan block with its dimension > 1 must be less than 1. 

Eq. (3.27) and its proof strongly suggest that the stability boundary of the a(3) scheme is defined by 
|^| = 1/2. In fact this conjecture has been verified numerically. We conclude this section by studying this 
special case. 

3.4. The |^| = 1/2 case 
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Let v = 1/2. Then Eq. (3.17) reduces to 


a 3 - e w a 2 - e~ i8 a + 1 = (a — e ie ) (a - e~ i8/2 ^ (a + e~ i8/2 ^j = 0 (^=1/2) (3.61) 

Thus the eigenvalues er are 

a = a±(9) d = ±e~ t0 ^ 2 and a = a 0 (9) d = e lS (u = 1/2) (3.62) 

On the other hand, by using Eqs. (3.14) and (3.62), one concludes that the eigenvalues a for the case 
v = —1/2 are 

a = a±(9) = ±e 10 / 2 and cr = a 0 (9) = e~ l8 ( v = —\j2 ) (3.63) 

For each of the above two cases, Eqs. (3.62) and (3.63) imply that the three eigenvalues are distinct if 

n = 0, ±1, ±2, . . . (3.64) 




Also, because the analytical amplification factor is e~ lv 8 for any (i/, 9) (see p.4 of [61]), for the case v = 1/2 
(v = —1/2), one of the roots of Eq. (3.17), i.e., cr+{9) (a + (9)), is identical to the analytical amplification 
factor. 

Consider the plane wave solution 

u(x,t) = e ik ^ x - at) 


The period associated with this solution is 


T = 


27T 

|fca| 


Let n, the number of total marching steps, and At be chosen such that 

nAt = NT, N = 1,2,3,... 

Then one has 

2irN 


\o\W\ 


where 


9 = kAx 

is the variation of phase angle over the interval ax. For the case |z^| = 1/2, Eq. (3.68) reduces to 

47 riV 


n = 


\ 0 \ 


(3.65) 

(3.66) 

(3.67) 

(3.68) 

(3.69) 

(3.70) 


Eqs. (3.62), (3.63), and (3.70) imply that 


[a ± {9)] n = a ± (0) =(±ir 


and 


Thus 


= 1 


[<r±m n = <T±W) 


koW] = I O'o(6 ) )j 

= 1, if n is even 


(3.71) 

(3.72) 

(3.73) 
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On the other hand, Eq. (3.68) implies that 


(e~ ive ) n = 1 (3.74) 

Eqs. (3.72)-(3.74) and the fact that e~ lv8 is the analytical amplification factor imply that the numerical 
solution generated by the a(3) scheme in a simulation involving a periodic boundary condition, aside from 
round-off errors, should be identical to the exact solution if n and At are chosen according to Eq. (3.67), 
and n is even, and (ii) the phase angles of the Fourier components involved in the simulation observe the 
condition Eq. (3.64) (see Sec. 5 in [1]). This prediction has been verified numerically (see Sec. 4). 

4. Numerical results 

To assess the accuracy of the a(3) scheme, consider the model problem with the PDE 


du 

at 



and the exact solution 


u = u e (x, t) = { sin [27r(a: — i)] = ^- T ^ 


gi27r(:r— t) ^—i27c(x—at) 


We have 


a=L = T= 1 


where L = wavelength and T = period. Let (i) 


, . def du e (x,t) . . def d 2 u e (x,t) 

u xe {x,t) = — — and u xxe (x,t) = 


dx 


dx 2 


(4.1) 


(4.2) 

(4.3) 


(4.4) 


and (ii) the spatial domain of unit length be divided into K uniform intervals. Thus 


ax = 1 /K, At = vax and t = nAt 


(4.5) 


where n = number of time steps, and t = total marching time. 

It has been shown numerically that the a(3) scheme is stable if 

M < 1/2 (4.6) 

On the other hand, (i) the a scheme is stable if 

M < 1 (4-7) 

and (ii) the a(4) scheme [72], another high order solver of Eq. (1.1), is stable if 

M < 1/3 (4.8) 

In Tables 1-4, the numerical errors of several computations using the a(3) and a schemes are presented 
in terms of the following error norms for the non-normalized mesh variables: 


E(K, n, v) 


def 


N 


l 

K 


K - 1 


Eb" 

3 = 0 


U e (Xj,t n )\ 2 


(4.9) 
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(4.10) 


and 


E x (K,n,u) = f 


\ 


1 K _1 

~ u xe{Xj,t n )] 2 

J=0 


E xx (K,n,u) d = 


N 


i * _i 

^ [(Uxx)" - Wxxe(a;j,i ")] 2 
A 3=0 


(4.11) 


The numerical errors of several simulations with v = 0.1 and t = 9.876 are given in Table 1. For the 
a scheme, as the values of K and n become larger, the values of E and E x are both reduced by a factor of 
about 4 as both K and n double their values, i.e., the scheme is 2nd order in accuracy for both it” and (w x )”. 
On the other hand, for the a(3) scheme, the values of E, E x , and E xx are reduced by the factors 16, 16, and 
4, respectively. Thus the a(3) scheme is 4th order in accuracy for both it” and ( u x )" while only 2nd order in 
accuracy for (tt xx )”. From the results shown, one can see that the a(3) scheme is much more accurate than 
the a scheme. As an example, for the case with K = 25 and n = 2469, the value of E for the a scheme is 
larger than that for the a(3) scheme by a factor of 3450! Because the a scheme is only 2nd order in accuracy 
for it", it is estimated that the accuracy of it” achieved by the a(3) scheme with K = 25 and n = 2469 is 
identical to that achieved by the a scheme with K = 25 x v / 3450 ss 1468 and n = 2469 x -\/3450 ss 145029. 

In Table 2, the cases considered have v = 0.1 and t — 10.00 = 10T. For these cases where t, is an integer 
multiple of the period T, it is seen that the a(3) scheme is 4th order in accuracy for it”, (w x )", and (it xx )”. 

In Table 3, the cases considered have v = 0.5 and t = 49.38, For these cases where the value of v is 
right at the stability boundary of the a(3) scheme, aside from round-off errors, the numerical values of (u x )” 
generated using the a(3) scheme are all identical to their exact solution values, respectively, if n are even 
integers. 

In Table 4, the cases considered have v = 0.5 and t = 50.00 = 50T, For these cases where (i) the 
value of v is right at the stability boundary of the a(3) scheme and (ii) t, is an integer multiple of T, aside 
from round-off errors, the numerical values of it”, (u x )”, and (it xx )" generated using the a(3) scheme are 
all identical to their exact solution values, respectively. Note that: (i) n and At are chosen according to 
Eq. (3.67) and n is even for each of these cases, and (ii) the exact solution are a linear combination of two 
plane wave solutions with \9\ = 27rAa: < 27r/25, i.e., 9 observes the condition Eq. (3.64). Thus the results 
shown in Table 4 confirm the prediction made at the end of Sec. 3. 


5. Conclusions and discussions 


A thorough and rigorous discussion of a new high order neutrally stable CESE solver of Eq. (1.1) has 
been presented. Because this two-level explicit scheme is associated with three independent mesh variables 
and three equations per mesh point, it is referred to as the a(3) scheme. As in the case of other similar CESE 
neutrally stable solvers [1,5,11,72], the a(3) scheme enforces conservation laws locally and globally, and it 
has the basic, forward marching, and backward marching forms. These forms are equivalent and satisfy the 
STI invariant property defined in Sec. 2. 

Based on the concept of STI invariance, a set of algebraic relations (Eqs. (2.114)-(2.118)) involving the 
coefficient matrices Qo(v), Q+(v) and Q-(v) is developed in Sec. 2. As it turns out, these relations can be 
used to construct a simple proof for the fact that the a(3) scheme is neutrally stable (i.e., non-dissipative) 
when it is stable. Another set of algebraic relations (Eq. (2.119)) involving the same matrices and its relation 
to other invariance property are also discussed in Sec. 2. 

In addition to establishing the neutral stability of a(3) scheme, in Sec. 3, it is also proved analytically 
that all three amplification factors of the a(3) scheme are of unit magnitude for all phase angles if \v\ < 1/2. 
This amazing theoretical result is completely consistent with the numerical stability condition \v\ < 1/2. 

It is shown in Sec. 4 that the a(3) scheme generally is (i) 4th-order accurate for the mesh variables it” 
and (it x )”; and 2nd-order accurate for (it xx )”. However, in some exceptional cases, the scheme can achieve 
perfect accuracy aside from round-off errors. The stability bound \is\ < 1/2 of the a(3) scheme is lower than 
the stability bound \v\ < 1 of the a scheme, but high than the stability bound \v\ < 1/3 of the a(4) scheme. 
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The CESE development has been driven by a basic idea that each practical scheme be built from a non- 
dissipative core scheme so that the numerical dissipation can be controlled effectively. As such, development 
of the a( 3) and a(4) schemes provides a solid foundation for the development of other more practical high 
order CESE schemes. 


Appendix A. Proof for Theorems 1 and 2 in Sec. 3 

First we prove Theorem 1. According to Eq. (8.20) on p.265 of [75], the fact that Xe, £ = 1,2 ,N 
are the eigenvalues of the N x N matrix A <t=> 

det(A - XI) = (Ai - A)(A 2 - A) • • • (X N - X) (A. 1) 

where A is any complex variable and / is the N x N identity matrix. Let A = 0. Eq. (A.l) implies that 


det (A) = AiA 2 • • • Ajv {A. 2) 

By definition, A is nonsingular <t=t det(A) y^ 0. Thus part (i) follows from Eq. (A. 2). 

According to Eq. (A.l), to prove part (ii) we need only to show that 

(£-*)•• ■(£■*) (A3) 

for any complex variable A. Because det(BC') = clet(13) det(C) for any two N x N matrices B and C, we 
have det(A)det (A -1 ) = det (AA _1 ) = det (/) = 1, i.e., det(A _1 ) = l/det(A). Thus Eq. (A. 2) implies that 


det (A 1 ) 


1 

A1A2 • • • Xn 


(A.4) 


By comparing Eqs. (A. 3) and (A.4), one concludes that Eq. (A. 3) is valid if A = 0. 
Let A y^ 0. Then 


A -1 - XI = -XI A- 1 



(A. 5) 


Thus 

A-il) (A.6) 

With the aid of (i) Eqs. (A.l) and (A.4) and (ii) the fact that det(— XI) = (—X) N , Eq. (A.6) implies that 


det (A 1 — A I) = det(— XI) det (A 1 ) det ^ 


det (A 1 — A I) 


(~X) N 










(A. 7) 


i.e., Eq. (A. 3) is also valid if A 0. Thus part (ii) of Theorem 1 has been proved. QED. 
According to Eq. (A.l), to prove Theorem 2, we need only to show that 


det (A - XI) = (Ai — A)(A 2 — A) • • • ( X N - X) 


(A. 8) 


Because det(M) = det(M), by using Eq. (A.l), we have 

det(A — XI) = det^A — A I^j = det(A — XI) 

= (Ai - A)(A 2 - A) • • • (Ajv - A) = (AT - A)(A^ - A) • • • (Xn - X) 


(A. 9) 
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i.e., Eq. (A. 8) has been proved. QED. 
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TABLE 1.— NUMERICAL RESULTS OF THE a( 3) AND a SCHEMES 


v = 0.1 

1 = 9.876 


77 = 25, n = 2,469 

77 =50, n = 4,938 

77= 100, n = 9,876 

77 =200, n= 19,752 

F 

fl( 3) 

0.131xl0“ 3 

0.143xl0“ 4 

0.883xl0“ s 

0.549xl0“ 7 


a 

0.452 

0.115 

0.287x10“' 

0.716xl0“ 2 

f 

«(3) 

0.445x10”’ 

0.977xl0“ 3 

0.611x10"* 

0.382xl0” 5 

C'X 

a 

2.90 

0.732 

0.182 

0.454x10“’ 

f 

C'xx 

0(3) 

0.225 

0.169 

0.406x10“’ 

0.100x10“’ 


TABLE 2.— NUMERICAL RESULTS OF THE a( 3) AND a SCHEMES 


v = 0.1 

1= 10.00 


77 = 25, n = 2,500 

77=50, n = 5,000 

77= 100, n = 10,000 

77= 200, n = 20,000 

F 

0(3) 

0.228xl0“ 3 

0.1 10xl0' 4 

0.628xl0“ s 

0.384xl0“ 7 


a 

0.469 

0.118 

0.292x10“’ 

0.727x1 0“ 2 

F 

0(3) 

0.154x10“’ 

0.992xl0“ 3 

0.623X10" 4 

0.390xl0“ 5 

C'x 

a 

2.89 

0.728 

0.182 

0.455x10“’ 

F 

^xx 

0(3) 

0.473 

0.316x10“’ 

0.199xl0“ 2 

0.124xl0“ 3 


TABLE 3.— NUMERICAL RESULTS OF THE a( 3) AND a SCHEMES 


v = 0.5 

7 = 49.38 


77 = 25, « = 2,469 

77=50, n = 4,938 

77= 100, « = 9,876 

77= 200, h= 19,752 

F 

0(3) 

0.168xl0“ 3 

0.471xl0“ 5 

0.294xl0“ 6 

0.183xl0“ 7 


a 

1.34 

0.429 

0.109 

0.271x10“’ 

F 

0(3) 

0.583x10“’ 

0.856xl0“ 12 

0.261x10“” 

0.678x10“” 


a 

8.73 

2.73 

0.686 

0.171 

F 

^xx 

0(3) 

1.01 

0.942x10’ 1 

0.235x10“’ 

0.587xl0“ 2 


TABLE 4.— NUMERICAL RESULTS OF THE a( 3) AND a SCHEMES 


v = 0.5 

7 = 50.00 


K = 25, n = 2,500 

77=50, « = 5,000 

77= 100, n = 10,000 

77= 200, n = 20,000 

F 

0(3) 

0.362xl0“ 13 

0.140xl0“ 12 

0.229xl0“ 12 

0.262xl0“ 12 


a 

1.35 

0.440 

0.111 

0.275x10“’ 

E 

0(3) 

0.162xl0“’ 2 

0.845 xl0“‘ 2 

0.261x10“” 

0.682x10“” 


a 

8.73 

2.73 

0.689 

0.172 

F 

^xx 

0(3) 

0.172xl0“ 9 

0.282xl0“ 8 

0.185xl0“ 7 

0.840xl0“ 7 
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